1 Sissejuhatus R keelde

1.1 Teaser

1.1.1 Paketid

library(tidyverse)

1.2 Taustaks

R on statistikale ja andmete visualiseerimisele suunatud vabavaraline programmeerimiskeel. R loodi 90ndate alguses S-i baasile ning on tasapisi muutunud üheks enimkasutatavaks statistika tarkvaraks.

R Passes SPSS in Scholarly Use, Stata Growing Rapidly R Passes SPSS in Scholarly Use, Stata Growing Rapidly (2014)

R-il on olulisi eeliseid mitmete teiste statistikatarkvarade ees:

  • R on eeskätt programmeerimiskeel, mis tähendab, et see võimaldab luua uusi funktsioone vastavalt eesseisvale probleemile;
  • R on vabavaraline ning kasvava populaarsusega, mis tähendab, et uute meetodid jõuavad R-i kasutajateni enne teisi;
  • R jookseb kõikidel levinud OS-idel ning on tasuta kättesaadav;
  • R ühildub mitmete teiste keeltega, mis teeb selle laiemalt kasutatavaks (näiteks käesolev markdown-i dokument);
  • R pakub visualiseerimisel tõenäoliselt parimat tasakaalu lihtsuse ja esteetilisuse vahel;
  • R on tõenäoliselt üks parimaid vahendeid statistika tegemiseks ka aastakümne pärast, mis tähendab, et selle õppimiseks kulutatud aeg pole raisatud.

Sellega seoses on sel ka miinuseid:

  • nõuab uue keele õppimist, mis võib olla alguses keeruline;
  • aeg ajalt esineb uuemates pakettides vigu, mis tähendab, et tuleb olla tähelepanelik;
  • seab kõrgemad nõuded meetodite tundmisele.

1.4 Installeerimiseks

R-i kasutamiseks piisab vaid R-i “mootori” installeerimisest, kuid mugavam on kasutada graafilist kasutajaliidest, näiteks R-Studiot.

Kuna R on vabavaraline tarkvara, siis on mitmete R-i funktsioonide kasutamiseks vajalik need esmalt alla laadida. See toimub R-i keskkonnas install.packages(...) abil (nt install.packages("tidyverse")), mispuhul pakett on funktsioonide kogum. Vajaminevate pakettide peale hetkel mõtlema ei pea, need tulevad jooksvalt praktikumide käigus.

1.5 Kui satud hätta..

R-i pakettidel on enamasti küllatki hästi lahti kirjutatud help sektsioon, mis on leitav help tabi alt (enamasti paremal all nurgas). Help tabi otsingusse saab kirjutada otsitud funktsiooni nime, misjärel ta kuvab ta funktsiooni kirjelduse ning argumentide väärtused. Enamasti on lisatud mõned näited funktsiooni kasutamisvõimalustest. Võimalik on veel:

  • kirjutada konsooli ?funktsiooni nimi (nt ?read.csv()) või help(...), kuhu läheb argumendiks funktsiooni nimi;
  • selekteerida funktsioon ja vajutada klaviatuuril F1.

Lisaks R-is olemasolevale kirjeldusele on mitmed paketid avaldatud teadusajakirjades, eeskätt ajakirjas Journal of Statistical Software. Samuti on enamus tüüp-probleemid lahendatud Stack Overflow lehel. Mõlemal juhul peaks google viima sobiva lahenduseni.

1.6 Tööpõhimõtted

1.6.1 Objektid

R-i põhiline töövahend on objekt, mis võib sisaldada nii funktsioone kui andmeid. Objekte saab ise luua <- või = märkide abil (<- on mõnevõrra kindlam meetod ja mõned soovitavad kasutada seda. Üldjuhul, sh selle aine raames, piisab ka = kasutamisest).

x = 5
x
[1] 5

Kui soovime siduda rohkem andmeid ühe objektiga, on võimalik luua vektor (jada / ühedimensionaalne maatriks).

x <- c(5, 12.3, 3:10)
x
 [1]  5.0 12.3  3.0  4.0  5.0  6.0  7.0  8.0  9.0 10.0

Objekt võib saada sisendi ka teisest objektidest. Samuti võib teha objektidega tehteid.

y = x - 5
y
 [1]  0.0  7.3 -2.0 -1.0  0.0  1.0  2.0  3.0  4.0  5.0

Objektiks võib olla ka funktsioon.

lahuta.viis = function(x){
  tmp = x - 5
  return(tmp)
}

y = lahuta.viis(x)
y
 [1]  0.0  7.3 -2.0 -1.0  0.0  1.0  2.0  3.0  4.0  5.0

1.6.2 Põhilised andmetüübid

Numbrid:

x = 1
str(x)
 num 1

Tähemärgid (numbrid muudetakse tähemärkideks):

x = c("a", "b", 1, 2)
str(x)
 chr [1:4] "a" "b" "1" "2"

Faktorid (identsed elemendid saavad sama koodi):

x = as.factor(c("a", "b", "1", "2", "a", "2"))
str(x)
 Factor w/ 4 levels "1","2","a","b": 3 4 1 2 3 2

Kuupäevad:

x = as.Date("2017-09-01")
str(x)
 Date[1:1], format: "2017-09-01"

Kuupäevade puhul on oluline jälgida formaati. Vaikimisi on selleks yyyy-mm-dd, kuid kui sisse loetavad andmed on teises formaadis, tuleks see as.Date(...) käskluses täpsustada format või origin argumentidega. Formaatidest pikemalt: Date Values.

as.Date("01-09-2017", format="%d-%m-%Y")
[1] "2017-09-01"
as.Date(17410, origin="1970-01-01")
[1] "2017-09-01"

1.6.3 Andmetabelid

Seni vaadatud jadad on olnud kõik ühedimensionaalsed, ent andmed võivad olla ka mitmedimensionaalsed. Enimlevinud mitmedimensionaalsed formaadid on andmetabelid ja maatriksid, millest peamiselt keskendume esimesele. Maatriks on kahedimensionaalne vektor, mis on justkui mitu vektorit kokku liidetuna, kuid kõik vektorid peavad olema sama liiki andmetega.

matrix(1:9, nrow=3, byrow = T)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
matrix(letters[1:9], nrow=3, byrow = T)
     [,1] [,2] [,3]
[1,] "a"  "b"  "c" 
[2,] "d"  "e"  "f" 
[3,] "g"  "h"  "i"

Mitu samapikkust vektorit saab liita kokku üheks ka andmetabeliks ehk data frameiks (data.frame(...)). Data frame’i erinevus seisneb selles, et veerud võivad sisaldada eritüüpi andmeid (nt üks veerg on numbriline, teine sisaldab faktoreid jne).

x = c(1:4)
y = c("Audi", "DAF", "Ford", "GAZ")
z = x-5

df = data.frame(tunnus.1 = x,
                tunnus.2 = y,
                tunnus.3 = z)
df
  tunnus.1 tunnus.2 tunnus.3
1        1     Audi       -4
2        2      DAF       -3
3        3     Ford       -2
4        4      GAZ       -1

Andmetabelite tunnuseid (ehk veerge) saab kätte ka vektorina kasutades $ märki:

df$tunnus.2
[1] Audi DAF  Ford GAZ 
Levels: Audi DAF Ford GAZ

Võimalik on kasutada ka indeksi (järjekorra) numbrit kandilistes sulgudes. Andmetabelite puhul on esimene number rea number ning koma järel veeru number. Vektori puhul piisab vaid rea numbrist.

df[,2]
[1] Audi DAF  Ford GAZ 
Levels: Audi DAF Ford GAZ
df["tunnus.2"]
  tunnus.2
1     Audi
2      DAF
3     Ford
4      GAZ
df[1,]["tunnus.2"]
  tunnus.2
1     Audi
df[1,2]
[1] Audi
Levels: Audi DAF Ford GAZ
df$tunnus.2[2]
[1] DAF
Levels: Audi DAF Ford GAZ

1.6.4 Listid

Kõiki erinevaid objekte saab koondada listidesse (list(...)). List on hea koondamaks erinevate pikkustega või erineva iseloomuga, kuid samateemalisi andmeid. Tihti võib liste kohata erinevate funktsioonide väljundites.

ls = list(df)
ls[[2]] = c(1:10)
ls[[3]] = y[2:4]

ls
[[1]]
  tunnus.1 tunnus.2 tunnus.3
1        1     Audi       -4
2        2      DAF       -3
3        3     Ford       -2
4        4      GAZ       -1

[[2]]
 [1]  1  2  3  4  5  6  7  8  9 10

[[3]]
[1] "DAF"  "Ford" "GAZ"

1.6.5 Funktsioonid objektide kohta

Kõikide objektide puhul:

  • str(...) - objekti struktuur
  • typeof(...) - storage mode
  • is.numeric(...) - kas on…
  • is.character(...)
  • is.factor(...)
  • is.function(...)
  • is.logical(...)

…

Vektori puhul:

  • length(...) - vektori pikkus

Data frame’i puhul:

  • dim(...) - data frame’i dimensioonid (ridu, veerge)
  • nrow(...) - ridade arv
  • ncol(...) - veergude arv

Formaadi muutmiseks:

  • as.vector(...) - vektor
  • as.matrix(...) - maatriks
  • as.Date(...) - kuupäev
  • as.data.frame(...) - data frame
  • as.list(...) - list

1.7 Andmestike laadimine

R suudab lugeda sisse enamikke levinud andmefailide formaate, kuid kõige sagedasemalt kasutatakse .csv (comma separated values). Andmestike sisse lugemiseks on erinevaid funktsioone, kuid .csv-de puhul piisab read.csv(...). read.csv(...) võimaldab muu hulgas määratleda, kuidas on eraldatud erinevad veerud (sep; koma, semikooloni vm abil) ja mida kasutatakse komakohtade eraldamiseks (dec; punkt, koma vm). See võib olla oluline, kuna erinevad programmid erinevatel OS-idel võivad kasutada erinevaid standardeid.

1.7.1 Internetist

url = "https://martenveskimae.github.io/stat-mod/party_dataset.csv"
df = read.csv(url)

Andmestiku vaatamiseks võib küsida tervet andmestikku, aga võib ka küsida esimesi või viimaseid ridu. Selleks saab kasutada nii indeksi kui ka eraldi funktsioone head(...) ja tail(...).

df[1:6, 1:4]
  party support support_emor  last
1   ref    0.36           NA 0.278
2   ref    0.34           NA 0.278
3   ref    0.45           NA 0.278
4   ref    0.43           NA 0.278
5   ref    0.43           NA 0.278
6   ref    0.44           NA 0.278
head(df[,1:4])
  party support support_emor  last
1   ref    0.36           NA 0.278
2   ref    0.34           NA 0.278
3   ref    0.45           NA 0.278
4   ref    0.43           NA 0.278
5   ref    0.43           NA 0.278
6   ref    0.44           NA 0.278
tail(df[,1:4])
    party support support_emor  last
551   sde      NA         0.15 0.152
552   irl      NA         0.06 0.137
553  kesk      NA         0.23 0.248
554   ref      NA         0.26 0.277
555  vaba      NA         0.07 0.087
556  ekre      NA         0.16 0.081

Et näha, millised tunnused andmestikus on, võib küsida tunnuste pealkirju colnames(...) funktsiooiga.

colnames(df)
 [1] "party"                 "support"              
 [3] "support_emor"          "last"                 
 [5] "year"                  "month"                
 [7] "gov"                   "month_cycle"          
 [9] "cycle_constant"        "cycle"                
[11] "inflats_m_basevalue"   "inflats_m"            
[13] "unemp"                 "gdpgr"                
[15] "gdpgr_yearly"          "skp_yearly"           
[17] "support_lag1"          "bruto"                
[19] "bruto_yearly"          "inflats_prevyear_base"
[21] "alko_suits_base"       "alko_suits_prevyear"  
[23] "eluase_base"           "eluase_prevyear"      
[25] "irl_meedias"           "kesk_meedias"         
[27] "ekre_meedias"          "ref_meedias"          
[29] "sde_meedias"           "rohel_meedias"        
[31] "date"

1.7.2 Kohalikult kettalt

Kohalikult kettalt andmestike laadimiseks tuleb märkida ära andmestiku asukoht kettal. Selle lihtsustamiseks on otstarbekas märkida ära töökeskonna asukoht.

Olemasoleva asukoha nägemiseks on funktsioon getwd().

getwd()
[1] "/Users/martenveskimae/GitHub/stat-mod"

Kui tahta seda muuta, saab kasutada setwd(...) funktsiooni, kus argumendiks on soovitud uus asukoht.

1.8 Andmetega tutvumine

Kokkuvõtvaid andmed nii andmetabelite kui teiste andmeformaatide kohta näeb funktsiooniga summary(...).

summary(df)
  party        support       support_emor         last       
 ekre: 30   Min.   :0.040   Min.   :0.0500   Min.   :0.0810  
 irl :124   1st Qu.:0.130   1st Qu.:0.1000   1st Qu.:0.1520  
 kesk:124   Median :0.210   Median :0.1600   Median :0.2100  
 ref :124   Mean   :0.204   Mean   :0.1601   Mean   :0.2014  
 sde :124   3rd Qu.:0.260   3rd Qu.:0.2200   3rd Qu.:0.2610  
 vaba: 30   Max.   :0.450   Max.   :0.2800   Max.   :0.2900  
            NA's   :34      NA's   :418      NA's   :12      
      year          month             gov          month_cycle
 Min.   :2007   Min.   : 1.000   Min.   :0.0000   Min.   : 1  
 1st Qu.:2010   1st Qu.: 3.000   1st Qu.:0.0000   1st Qu.:11  
 Median :2012   Median : 6.000   Median :1.0000   Median :20  
 Mean   :2012   Mean   : 6.378   Mean   :0.5665   Mean   :22  
 3rd Qu.:2015   3rd Qu.: 9.000   3rd Qu.:1.0000   3rd Qu.:33  
 Max.   :2017   Max.   :12.000   Max.   :1.0000   Max.   :49  
                                                  NA's   :12  
 cycle_constant     cycle         inflats_m_basevalue   inflats_m       
 Min.   :49     Min.   :0.02041   Min.   :108.5       Min.   :-0.90701  
 1st Qu.:49     1st Qu.:0.22449   1st Qu.:124.7       1st Qu.:-0.03494  
 Median :49     Median :0.40816   Median :138.0       Median : 0.28775  
 Mean   :49     Mean   :0.44898   Mean   :134.2       Mean   : 0.30774  
 3rd Qu.:49     3rd Qu.:0.67347   3rd Qu.:143.9       3rd Qu.: 0.64917  
 Max.   :49     Max.   :1.00000   Max.   :146.3       Max.   : 2.08868  
 NA's   :12     NA's   :12        NA's   :72          NA's   :180       
     unemp            gdpgr            gdpgr_yearly     
 Min.   : 4.000   Min.   :-10.80000   Min.   :-18.4000  
 1st Qu.: 6.500   1st Qu.: -0.30000   1st Qu.: -1.7000  
 Median : 8.000   Median :  0.70000   Median :  2.3000  
 Mean   : 9.141   Mean   : -0.01209   Mean   :  0.7391  
 3rd Qu.:11.300   3rd Qu.:  1.30000   3rd Qu.:  5.6000  
 Max.   :19.500   Max.   :  3.70000   Max.   : 11.6000  
 NA's   :78       NA's   :192         NA's   :162       
   skp_yearly        support_lag1        bruto         bruto_yearly   
 Min.   :-19.3000   Min.   :0.0700   Min.   : 659.7   Min.   :-6.500  
 1st Qu.:  0.3000   1st Qu.:0.1500   1st Qu.: 792.3   1st Qu.: 3.900  
 Median :  1.9000   Median :0.2300   Median : 857.1   Median : 6.300  
 Mean   :  0.8247   Mean   :0.2231   Mean   : 898.4   Mean   : 6.707  
 3rd Qu.:  5.2000   3rd Qu.:0.2800   3rd Qu.:1010.0   3rd Qu.: 8.500  
 Max.   : 10.5000   Max.   :0.4500   Max.   :1163.0   Max.   :21.166  
 NA's   :78         NA's   :184      NA's   :84       NA's   :192     
 inflats_prevyear_base alko_suits_base alko_suits_prevyear  eluase_base   
 Min.   :102.7         Min.   :106.8   Min.   :102.1       Min.   :121.7  
 1st Qu.:123.2         1st Qu.:138.1   1st Qu.:115.3       1st Qu.:149.6  
 Median :133.1         Median :149.3   Median :141.2       Median :168.3  
 Mean   :130.6         Mean   :148.9   Mean   :139.5       Mean   :168.3  
 3rd Qu.:143.3         3rd Qu.:167.3   3rd Qu.:158.3       3rd Qu.:191.4  
 Max.   :145.7         Max.   :182.4   Max.   :174.9       Max.   :199.0  
 NA's   :60            NA's   :174     NA's   :174         NA's   :126    
 eluase_prevyear  irl_meedias      kesk_meedias     ekre_meedias  
 Min.   :106.1   Min.   : 227.0   Min.   : 312.0   Min.   :  8.0  
 1st Qu.:144.7   1st Qu.: 438.5   1st Qu.: 559.8   1st Qu.: 29.0  
 Median :158.0   Median : 559.0   Median : 626.5   Median : 51.5  
 Mean   :159.5   Mean   : 610.4   Mean   : 717.0   Mean   :115.7  
 3rd Qu.:183.2   3rd Qu.: 730.2   3rd Qu.: 811.5   3rd Qu.:111.0  
 Max.   :199.0   Max.   :1597.0   Max.   :2414.0   Max.   :673.0  
 NA's   :126     NA's   :316      NA's   :316      NA's   :316    
  ref_meedias      sde_meedias     rohel_meedias            date    
 Min.   : 343.0   Min.   : 245.0   Min.   :  8.00   2015-01-01:  6  
 1st Qu.: 537.8   1st Qu.: 367.0   1st Qu.: 26.75   2015-02-01:  6  
 Median : 658.5   Median : 444.5   Median : 49.00   2015-03-01:  6  
 Mean   : 747.7   Mean   : 510.6   Mean   : 85.40   2015-04-01:  6  
 3rd Qu.: 929.0   3rd Qu.: 585.8   3rd Qu.:103.00   2015-05-01:  6  
 Max.   :1836.0   Max.   :1480.0   Max.   :462.00   2015-06-01:  6  
 NA's   :316      NA's   :316      NA's   :316      (Other)   :520

Faktortunnuste puhul võib pakkuda huvi, millised erinevad kategooriad on andmestikus esindatud. Näiteks võime vaadata, millised erakonnad andmestikus leiduvad.

levels(df$party)
[1] "ekre" "irl"  "kesk" "ref"  "sde"  "vaba"

Kui soovida täpsemat ülevaadet, näiteks, mitu andmepunkti on iga erakonnaga seotud, siis selleks on hea kasutada paketti tidyverse. Tidyverse hõlmab endas mitmeid erinevaid, kuid omavahel ühilduvaid pakette. Üks neist, dplyr, pakub mitmeid funktsioone andmete mugavaks töötlemiseks (eeskätt filtreerimine, grupeerimine, tunnuste muutmine jne). Teine oluline eelis tuleneb paketist magrittr, mis võimaldab käsklusi sisestada jadana (ing.k pipe). See kiirendab tööprotsessi ja muudab R-i koodi kergemini jälgitavaks. Pipe’ide kasutamiseks tuleb kahe funktsiooni vahel lisada %>% sümbolid. Enne näite juurde liikumist on vajalik tidyverse aktiveerida, milleks on funktsioon library(...).

library(tidyverse)

Pipe’i näide võib olla järgnev:

c(1:10) %>% lahuta.viis()
 [1] -4 -3 -2 -1  0  1  2  3  4  5

Järgnevalt tutvume, mitu rida on seotud üksikute erakondadega (sagedustabel).

df %>%
  group_by(party) %>%
  summarize(n = n())
# A tibble: 6 x 2
   party     n
  <fctr> <int>
1   ekre    30
2    irl   124
3   kesk   124
4    ref   124
5    sde   124
6   vaba    30

dplyr kohta leiab rohkem informatsiooni siit: Grammar of Data Manipulation • dplyr

1.9 Filtreerimine & joonised

Kõige kiirem viis andmete visualiseerimiseks on kasutada funktsioone plot(...) või hist(...), kus tuleb määratleda telje/telgede väärtused.

plot(df$year, df$support)

hist(df$support)

Paraku ei ole antud funktsioonid väga mugavad tidyverse’is kasutamiseks ning kujunduse muutmine on piiratud võimalustega. Kõige enim levinud viis andmete visualiseerimiseks R-i keskkonnas on paketiga ggplot2, mis on samuti tidyverse’i osa.

df %>%
  filter(party=="ref") %>%
  ggplot() +
  geom_point(aes(year, support))

df %>%
  filter(party=="ref") %>%
  ggplot() +
  geom_histogram(aes(support))

Võrdlusoperaatorid

Filtreerimisel või muude võrdluste tegemisel kasutatakse teatud tingmärke, mida nimetatakse võrldusoperaatoriteks. R-is on need järnevad:

Operaator Tulemus
== Võrdne
!= Mittevõrdne
> Suurem
< Väiksem
>= Suurem-võrdne
<= Väiksem-võrdne
%in% Sisaldab

Loogilised operaatorid

Olukordades, kus soovime mitut võrdlust/tingimust kombineerida, tuleb erinevad komponendid siduda omavahel loogiliste operaatoritega:

Operaator Tulemus
& Loogiline JA üksikute vektori komponentide kohta (väljund on TRUE/FALSE kõigi vektori komponentide kohta)
&& Loogiline JA kõigi vektori komponentide kohta (väljund on ainult üks TRUE/FALSE kogu vektori kohta)
| Loogiline VÕI üksikute vektori komponentide kohta (väljund on TRUE/FALSE kõigi vektori komponentide kohta)
|| Loogiline VÕI kõigi vektori komponentide kohta (väljund on ainult üks TRUE/FALSE kogu vektori kohta)
! Loogiline mittesamaväärsus

Näiteid: R - Operators

Kujunduse kiireks muutmiseks on erinevad valmislahendused, algusega theme_.

df %>%
  filter(party=="ref") %>%
  ggplot() +
  geom_point(aes(year, support)) +
  theme_bw()

ggplot2 kohta leiab rohkem informatsiooni siit: Elegant Data Visualisations Using the Grammar of Graphics

2 Jaotus ja keskmised

2.0.1 Paketid

library(tidyverse)
library(e1071)
library(reshape2)

2.0.2 Andmestik

df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")

2.1 Histogramm

Vaatame jaotuste järsakust ja sümmeetrilisust visuaalselt histogrammi abil. geom_histogram(...) annab vaikimisi sagedusjaotuse, kuid on võimalik valida ka protsente või tõenäosustihedust.

df %>%
  ggplot() +
  geom_histogram(aes(support, ..count../sum(..count..))) +
  geom_line(aes(support, ..density../100), stat = "density") +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  labs(y="protsent") +
  theme_bw()

ggplot võimaldab luua eraldi joonised erinevate gruppide jaoks kasutades facet_wrap(...) käsklust.

df %>%
  ggplot() +
  geom_histogram(aes(support, ..count../sum(..count..))) +
  geom_line(aes(support, ..density../100), stat = "density") +
  facet_wrap(~party) +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  labs(y="protsent") +
  theme_bw()

Silmas tuleb pidada, et sel puhul kalkuleerib ggplot tulpade laiused ja hulga üle kõikide kuvatud jooniste, mistõttu erinevad need üksikult kuvatud joonistest. Tulpade hulga ja laiuste kontrollimiseks on argumendid bins ja binwidth, millest esimene määrab tulpade hulga ja teine laiused. Neid käsklusi ei saa korraga kasutada (binwidth kirjutab binsi üle).

df %>%
  ggplot() +
  geom_histogram(aes(support, ..count../sum(..count..)), bins=10, binwidth=0.03) +
  geom_line(aes(support, ..density../100), stat = "density") +
  facet_wrap(~party) +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  labs(y="protsent") +
  theme_bw()

2.2 Järsakuse- ja assümmeetriakordaja

Vaatame järsakust ja sümmeetrilisust ka numbriliselt nn järsakuse- ja assümmeetriakordaja abil. Seda võimaldab pakett e1071, kust leiame käsklused skewness(...) ja kurtosis(...) (rakendatavad valemid on leitavad helpist). Kuna soovime rakendada antud käsklusi kõigi gruppide (erakondade) jaoks eraldi, siis tuleks andmestik esmalt filtreerida grupeerimistunnuse alusel (party) ning alles seejärel rakendada antud käsklusi.

df %>%
  group_by(party) %>%
  summarise(skewness = skewness(support, na.rm=T),
            kurtosis = kurtosis(support, na.rm=T))
# A tibble: 6 x 3
   party    skewness   kurtosis
  <fctr>       <dbl>      <dbl>
1   ekre -0.23712096 -0.8569726
2    irl -0.61654227 -0.2531924
3   kesk  0.29345645 -0.2409381
4    ref  0.30822221 -0.5851121
5    sde  0.44024113 -1.2410302
6   vaba  0.01834276 -1.3350304

Alternatiiv dplyr meetodile oleks for loop. Loopimine on küll võrdlemisi lihtne, kuid vahel mitte kõige efektiivsem viis:

parties = as.character(unique(df$party))
s = c()
k = c()
for(i in 1:length(parties)){
  s[length(s)+1] = skewness(df$support[df$party==parties[i]], na.rm=T)
  k[length(k)+1] = kurtosis(df$support[df$party==parties[i]], na.rm=T)
}
names(s) = parties
names(k) = parties
s;k
        ref        kesk         irl         sde        vaba        ekre 
 0.30822221  0.29345645 -0.61654227  0.44024113  0.01834276 -0.23712096 
       ref       kesk        irl        sde       vaba       ekre 
-0.5851121 -0.2409381 -0.2531924 -1.2410302 -1.3350304 -0.8569726

Sama tulemuse saab efektiivsemalt kasutades apply(...) käsklust (antud juhul sapply(...). Rohkem applyst: Tutorial on the R Apply Family), mis jäljendab eelneva meetodi loogikat.

s = sapply(as.character(unique(df$party)), function(x) skewness(df$support[df$party==x], na.rm=T))
k = sapply(as.character(unique(df$party)), function(x) kurtosis(df$support[df$party==x], na.rm=T))
s;k
        ref        kesk         irl         sde        vaba        ekre 
 0.30822221  0.29345645 -0.61654227  0.44024113  0.01834276 -0.23712096 
       ref       kesk        irl        sde       vaba       ekre 
-0.5851121 -0.2409381 -0.2531924 -1.2410302 -1.3350304 -0.8569726

2.3 Keskmine, mediaan ja mood

Moodi, ehk kõige sagedamini esineva väärtuse leidmiseks eraldi käsklus puudub, ent seda on lihtne leida sagedustabelist.

df %>%
  group_by(party, support) %>%
  summarise(n = n()) %>%
  na.omit() %>%
  group_by(party) %>%
  summarise(mode = support[n==max(n)][1])
# A tibble: 6 x 2
   party  mode
  <fctr> <dbl>
1   ekre  0.10
2    irl  0.16
3   kesk  0.26
4    ref  0.25
5    sde  0.12
6   vaba  0.11

dplyri saab kasutada ka keskmise ja mediaani leidmiseks.

df %>%
  group_by(party) %>%
  summarise(mean = mean(support,na.rm=T),
            median = median(support,na.rm=T))
# A tibble: 6 x 3
   party       mean median
  <fctr>      <dbl>  <dbl>
1   ekre 0.10952381   0.10
2    irl 0.13875000   0.15
3   kesk 0.25858333   0.26
4    ref 0.28716667   0.29
5    sde 0.16633333   0.14
6   vaba 0.09904762   0.10

Viimaks võib kanda leitud väärtused histogrammile, et neid visuaalselt kõrvutada. Kuigi ggplot oskab lisada joonistele funktsioone, ei suuda ta praeguses versioonis kuvada eraldi funktsioone eraldi gruppidele (nt kasutades facet_wrap(~party)). Otstarbekas on sellisel juhul teha uus andmetabel keskmiste väärtustega. Andmetabeli mugavamaks kasutamiseks pöörame tabeli laiast formaadist pikka formaati, kasutades selleks reshape2 paketti ja melt(...) funktsiooni. Rohkem reshape2 ja laia/pika formaadi kohta: An Introduction to reshape2.

keskmised = df %>%
  group_by(party, support) %>%
  summarise(n = n()) %>%
  na.omit() %>%
  group_by(party) %>%
  summarise(keskmine = mean(support),
            mediaan = median(support),
            mood = support[n==max(n)][1]) %>%
  melt("party",c("keskmine", "mediaan", "mood"),"Keskmised")

df %>%
  ggplot() +
  geom_histogram(aes(support)) +
  geom_line(aes(support, ..density..), stat = "density") +
  geom_vline(data=keskmised,aes(xintercept=value, linetype=Keskmised, color=Keskmised)) +
  facet_wrap(~party) +
  scale_x_continuous(labels=scales::percent) +
  theme_bw()

Lisaks mediaanile võivad huvipakkuvad olla ka teised protsentiilid (mediaan on 50 protsentiil), eeskätt 25 ja 75 (teisisõnu I ja III kvartiil). Veel võib aeg-ajalt kohta mõistet kvartiilhaare (interquartile range - IQR), mis tähistab vahemikku I ja III kvartiili vahel. R-is on protsentiilide leidmiseks käsklus quantile(...).

sapply(as.character(unique(df$party)), function(x) quantile(df$support[df$party==x],c(.025,.25,.5,.75,.975),na.rm=T))
          ref   kesk     irl     sde  vaba ekre
2.5%  0.17975 0.2000 0.04975 0.08000 0.065 0.07
25%   0.23000 0.2375 0.12000 0.11750 0.080 0.10
50%   0.29000 0.2600 0.15000 0.14000 0.100 0.10
75%   0.33000 0.2800 0.16000 0.23000 0.120 0.13
97.5% 0.43000 0.3300 0.20025 0.28025 0.135 0.14

Eelnevalt toodud protsentiilid saab graafida boxploti abil, mis kuvab veel lisaks erindid (outliers).

df %>%
  ggplot() +
  geom_boxplot(aes(party,support)) +
  scale_y_continuous(labels=scales::percent) +
  theme_bw()

3 Usaldusvahemikud

3.0.1 Paketid

library(tidyverse)
library(reshape2)

3.0.2 Andmestik

df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")

3.1 Standardhälve

Eelmisest korrast teame, et keskmine toetusprotsent erakonniti erineb, aga erineb ka kõikumine keskmise ümber. Osade erakondade puhul on keskmine toetus stabiilselt ühe tüüpväärtuse ümber, aga teistel mitte. Selle hindamiseks kasutatakse tüüpiliselt variatsiooni ja standardhälvet, mis on R-is leitav käsklustega var(...) ja sd(...).

Variatsioon on eeldatav erinevus tõelise väärtuse ja keskmise vahel ruudus:

\(Var(X) = E[(X-\mu)^2]\)

ehk teisisõnu erinevuste keskmine:

\(Var(X) = \frac{ \sum(X-\mu)^2}{N-1}\)

X = 1:10

var(X)
[1] 9.166667
sum((X - mean(X))^2) / (length(X)-1)
[1] 9.166667

(Miks on jagaja \(N-1\) ja mitte \(N\)? Vahe tuleneb sellest, kas meil on andmed kogu populatsiooni kohta või valim populatsioonist, mis on alati väiksem kui N!)

Variatsiooni leidmiseks on erinevused võetud ruutu. Et jõuda tagasi algskaala juurde, tuleks variatsioonist võtta ruutjuur, mis ongi standardhälve:

\(\sigma=\sqrt{E[(X-\mu)^2]}\)

X = 1:10

sqrt(var(X))
[1] 3.02765
sd(X)
[1] 3.02765

Standardhälvet võib seega mõista kui keskmist erinevust keskmisest.

Vaatame erakondade keskmisi ja standardhälvet lähemalt:

keskmised = df %>%
  group_by(party) %>%
  summarise(keskmine = mean(support, na.rm=T),
            s.ylemine = keskmine + sd(support, na.rm=T),
            s.alumine = keskmine - sd(support, na.rm=T)) %>%
  melt("party",c("keskmine", "s.ylemine", "s.alumine"), "Varieeruvus")

df %>%
  ggplot() +
  geom_line(aes(as.Date(date), support)) +
  geom_hline(data=keskmised,aes(yintercept=value, linetype=Varieeruvus, color=Varieeruvus)) +
  facet_wrap(~party) +
  scale_y_continuous(labels=scales::percent) +
  labs(x="") +
  theme_bw()

3.2 Standardviga ja usaldusvahemikud

Millest see kõikumine keskmine ümber tuleneb on meile teada. Erakondade toetus varieerub ajas erinevalt ja meie arvutatud keskmine näitab tüüpväärtust üle aja. Siin on aga veel üks aga! Erakonna igakuised toetusnumbrid on ise sammuti üks statistiline hinnang (enamasti saadud umbes 1000se valimiga küsitlusest)! Milline on erakondade tõeline toetus populatsioonis on kindlaks tehtav vaid kogupopulatsiooni küsitledes ja seda me teha ei saa. Sellest tulenevalt on meie küsitlustulemustes alati mingi viga sees. Mida saame teha, on hinnata toetusprotsendi vahemikku, kus tõeline toetusprotsent tõenäoliselt asetseks.

Teoreetiliselt tähendaks see olukorda, kus võtaksime populatsioonist, kus on tegelikud teotused erakondadele, lõpmatu arv juhuvalimeid ja leiaks nende valimite toetusprotsendid. Enamasti oleks tulemus sarnane tegeliku toetusega (mida me ei tea), samas teatud olukordades oleks tulemus ühele või teisele poole väga mööda. Kui me teeks nendelt juhuvalimitelt toetusprotsentide graafi, tuleks see normaaljaotusega graaf: 95% juhtumitest langeks +/-1.96 z-skoori vahele. Sellest tulenevalt saaksime hinnata, millisesse vahemikku jääb tõeline tulemus (95% tõenäosusega).

Praktikas seda teha ei saa, sest küsitlused on kallid. Küll aga on võimalik tuletada usaldusvahemik kasutades selleks standardviga. Standardviga kirjeldab, kui lähedal on valimi keskmine populatsiooni keskmisele, kasutades selleks üldjuhul variatsiooni ja valimi suurust. Standardvea leidmiseks on erinevaid meetodeid, sõltuvalt andmete iseloomust. Antud juhul on tegu proportsioonidega (erakonna toetust väljendatakse protsendina populatsioonist), mispuhul kasutame proportsiooni standardviga:

\(SE_p = \sqrt{\frac{p(1-p)}{N}}\)

kus:

  • p on protsent (antud juhul erakonna toetusprotsent)
  • N valimi suurus (antud juhul 1000)

Usaldusvahemiku leidmiseks tuleb korrutada standardviga soovitud z-skoori väärtusega (nt 95% usaldusvahemiku leidmiseks on see 1.96) ning liita või lahutada erakonna toetusprotsendist:

\(p \pm ci_{95} = p \pm 1.96 \cdot SE\)

Kuivõrd tegu on küllalt spetsiifilise küsimusega, siis on otstarbekas teha oma funktsioon.

ci = function(p,n=1000,z=1.96){
  list(u = p + (z*(sqrt(p*(1-p)/n))),
       l = p - (z*(sqrt(p*(1-p)/n))))
}

df %>%
  mutate(ylemine =  ci(support)[["u"]],
         alumine =  ci(support)[["l"]]) %>%
  melt(c("party", "date"), c("support", "ylemine", "alumine"),"Toetus") %>%
  ggplot() +
  geom_line(aes(as.Date(date), value, color=Toetus)) +
  facet_wrap(~party) +
  scale_y_continuous(labels=scales::percent) +
  labs(x="",y="protsent") +
  theme_bw()

Lihtsamaks lugemiseks võime jooned üle siluda.

df %>%
  mutate(ylemine =  ci(support)[["u"]],
         alumine =  ci(support)[["l"]]) %>%
  melt(c("party", "date"), c("support", "ylemine", "alumine"),"Toetus") %>%
  ggplot() +
  geom_smooth(aes(as.Date(date), value, color=Toetus), se=F, method="lm", formula= y~poly(x, 7), size=.7) +
  facet_wrap(~party) +
  scale_y_continuous(labels=scales::percent) +
  labs(x="",y="protsent") +
  theme_bw()

Alternatiivne lähenemine eelnevalt kirjeldatud meetodile oleks simuleerida erakondade toetusprotsente bootstrap meetodil. Bootstrapi idee seisneb olemasolevast valimist juhuslike valimite võtmises ning nende põhjal huvipakkuvate statistikute arvutamises (eeskätt variatsioon ja standardviga).

Kuna antud juhul on tegu proportsioonidega ning me soovime igat üksikut vaatlust eraldi simuleerida, siis võime seda teha mündiviske loogikal. Mündiviske puhul simuleerib R binaarset jaotust, mille keskmine on meie poolt antud ning jaotus on binominaalne. Korrates mündiviske katset piisava arvu kordi (antud juhul katsetame 500ga), saame arvutada simulatsiooni standardhälve ja standardvea.

Mündiviske simulatsiooniks kasutame paketti mosaic ning funktsioone do(...), rFlip(...) ja confint(...). Käesolevas katses: do(1000) kordab katset 1000 korda ning rFlip(party.n,x) viskab ühes katses münti nii mitu korda, kui palju on erakonnaga seotud vaatlusi, tõenäosusega x. confint(...) arvutab simulatsiooni usaldusvahemikud.

toetus.sim = function(x,party){
  require(mosaic)
  if(is.na(x)) return(NA)
  party.n = nrow(df[df$party==party & !is.na(df$support),])
  tmp = do(1000) * rflip(party.n,x)
  ci = confint(tmp$prop)
  list(u = ci[["97.5%"]],
       l = ci[["2.5%"]])
}

df %>%
  rowwise() %>%
  mutate(ylemine = toetus.sim(support,party)[["u"]],
         alumine = toetus.sim(support,party)[["l"]]) %>%
  melt(c("party", "date"), c("support", "ylemine", "alumine"),"Toetus") %>%
  ggplot() +
  geom_smooth(aes(as.Date(date), value, color=Toetus), se=F, method="lm", formula= y~poly(x, 7), size=.7) +
  facet_wrap(~party) +
  scale_y_continuous(labels=scales::percent) +
  labs(x="",y="protsent") +
  theme_bw()

4 Keskmiste võrdlus

4.0.1 Paketid

library(tidyverse)
library(effsize)

4.0.2 Andmestik

df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")
paneel = read.csv("https://martenveskimae.github.io/stat-mod/ep2014_paneel.csv")

4.1 T-test

Vaatame Eesti poliitika sukka ja saabast ehk Reformierakonda ja Keskerakonda eelmise valimistsükli ajal.

df %>%
  filter(party %in% c("ref", "kesk"),
         year>=2011 & year<2015) %>%
  ggplot() +
  geom_point(aes(as.Date(date),support)) +
  facet_wrap(~party) +
  scale_y_continuous(labels=scales::percent) +
  labs(x="") +
  theme_bw()

Nende keskmine toetus on antud perioodil:

df %>%
  filter(party %in% c("ref", "kesk"),
         year>=2011 & year<2015) %>%
  group_by(party) %>%
  summarise(keskmine = mean(support,na.rm=T))
# A tibble: 2 x 2
   party  keskmine
  <fctr>     <dbl>
1   kesk 0.2444681
2    ref 0.2772340

Tekib küsimus, kas Keskerakonna keskmine toetus üldse erineb oluliselt Reformi keskmisest toetusest? Kahe keskmise võrdlemiseks saab kasutada t-testi meetodit (R-is t.test(...)), mis hindab, kas kaks keskmist pärinevad samast populatsioonist või on tegu erinevate gruppidega. Lihtsustatult öeldes, võrreldakse t-testis keskmiste erinevust kogu variatsiooniga, ning hinnatakse saadud tulemuse tõenäosust (kas tegu oli tõenäolise või ebatõenäolise tulemusega, ehk statistiliselt oluline või mitte). T-testi eelduseks on, et mõlemad valimid on normaaljaotuslikud ning sarnase variatiivsusega (viimast eeldust on võimalik rikkuda kasutades “robustset” t-testi, kus vabadusastmeid arvutatakse teisel meetodil!; R kasutab vaikimisi robustset versiooni - argument var.equal = FALSE).

T-statistiku arvutamiseks ühe valimi puhul on valemiks:

\(t=\frac{\bar{x}-\mu_o}{\frac{s}{\sqrt{N}}}\)

kus

  • \(\mu_0\) on võrreldav keskmine
  • s on valimi standardhälve

..kahe sõltumatu valimi puhul:

\(t=\frac{\bar{x_1}-\bar{x_2}}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}}\)

..ja kahe sõltuva valimi puhul:

\(t=\frac{\frac{\sum{d}}{N}}{\sqrt{\frac{\sum{d^2}-\frac{(\sum{d^2})}{N}}{N(N-1)}}}\)

kus

  • d on paaris olevate väärtuste erinevus
  • N on paaride arv

Teeme meetodi illustreerimiseks simulatsiooni. Ütleme, et meil on kaks sõltumatut normaaljaotusega vektorit, mõlemad 50 vaatlusega. Arvutame vektorite t-testi tulemused ja kordame katset 1000 korda (iga kord uute andmetega).

t.statistikud = replicate(1000, t.test(rnorm(50), rnorm(50))$statistic)

Vaatame saadud t-statistikute jaotust visuaalselt

ggplot() +
  geom_density(aes(t.statistikud)) +
  theme_bw()

Kui meie eeldused peavad paika, siis peaks tulemus lähenema t-jaotusele vabadusastmel 98 (50+50-2). Genereerime t-jaotuse antud vabadusastmetega ja arvutame 95% usaldusvahemikud. T-jaotuse leidmiseks peaksime leidma t-statistikute umbkaudsed otspunktid, näiteks käsklusega range(...), seejärel looma sümmeetrilise (sest t-jaotus on sümmeetriline!) vektori käsklusega seq(...) ning arvutama saadud vektori t-jaotusliku tiheduse etteantud vabadusastmete puhul (praegu 98), kasutades käsklust dt(...). Nendest saavad teoreetilise jaotuse x ja y atribuudid. Usaldusvahemike leidmiseks kasutame käsklust qt(..), mis arvutab teoreetilise jaotuse kvantiilid (praegu 0.025 ja 0.975).

Lisame tulemused eelnenud graafile.

r = mean(abs(range(t.statistikud)))
x = seq(-r,r,length=100)
y = dt(x,df=98)

a95 = qt(c(.025, .975), df=98)

ggplot() +
  geom_density(aes(t.statistikud,color="simuleeritud t-statistikud")) +
  geom_line(aes(x,y,color="teoreetiline t-jaotus")) +
  geom_vline(xintercept = a95[1]) +
  geom_vline(xintercept = a95[2]) +
  theme_bw()

Joonte vahele jäänud tulemused olid statistiliselt ebaolulised ning servadesse sattunud tulemused statistiliselt olulised 95% usaldusnivoo puhul. Võime arvutada, mitu protsenti osutus olulisteks.

length(t.statistikud[t.statistikud<=a95[1] | t.statistikud>=a95[2]]) / length(t.statistikud)
[1] 0.043

4.1.1 Ühe valimi t-test

Hindame seda esmalt ühe valimi t testiga. Leiame Reformierakonna keskmise toetuse perioodil 2011-2014 ja võrdleme seda keskerakonna toetusega.

ref.keskmine = df %>%
  filter(party == "ref",
         year>=2011 & year<2015) %>%
  pull(support) %>%
  na.omit() %>%
  mean()

kesk = df %>%
  filter(party == "kesk",
         year>=2011 & year<2015) %>%
  pull(support) %>%
  na.omit()

t.test(kesk, mu=ref.keskmine)

    One Sample t-test

data:  kesk
t = -8.034, df = 46, p-value = 2.635e-10
alternative hypothesis: true mean is not equal to 0.277234
95 percent confidence interval:
 0.2362586 0.2526775
sample estimates:
mean of x 
0.2444681

4.1.2 Kahe sõltumatu valimi t-test

Kas valitsuse ja opositsiooni toetus erineb üldse oluliselt erinevatel ajahetkedel valimistsüklis? Selleks peaksime tegema kahe sõltumatu valimi t-testi. Testime seda esialgu vaid ühel aastal, näiteks 2014.

df %>%
  filter(year==2014) %>%
  select(support,gov) %>%
  t.test(support~gov,data=.)

    Welch Two Sample t-test

data:  support by gov
t = -2.3091, df = 45.844, p-value = 0.0255
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.055047311 -0.003770081
sample estimates:
mean in group 0 mean in group 1 
      0.2173913       0.2468000

Vaatame, kuidas t-testi p-väärtus varieerub erinevatel aastatel.

plyr::ldply(unique(df$year), function(x){
  data.frame(p.value = df %>%
               filter(year==x) %>%
               select(support,gov) %>%
               t.test(support~gov,data=.) %>%
               .$p.value,
             year = x)
  }) %>%
  ggplot() +
  geom_point(aes(year,p.value)) +
  geom_hline(yintercept = 0.05) +
  scale_x_continuous(breaks=2007:2017) +
  theme_bw()

Paneme lisaks kõrvale keskmised toetusnumbrid antud aastatel.

df %>%
  group_by(year,gov=as.factor(gov)) %>%
  summarise(support = mean(support,na.rm=T)) %>%
  ggplot() +
  geom_point(aes(year,support,color=gov)) +
  geom_line(aes(year,support,color=gov)) +
  scale_x_continuous(breaks=2007:2017) +
  scale_y_continuous(labels=scales::percent) +
  theme_bw()

4.1.3 Kahe sõltuva valimi t-test

Kahe sõltuva valimi t-testi erakonnatoetuse andmestikuga teha ei saa, sest meie uuringu disain ei vasta sellele. Selleks avame ühe teise andmestiku, milleks on paneeluuring enne ja pärast 2014.a. EP valimisi.

summary(paneel)
       t1.date    t1.evalimiste.usaldus       t2.date   
 2014-05-20:115   Min.   : 0.00         2014-06-05:222  
 2014-05-12:110   1st Qu.: 5.00         2014-06-03:183  
 2014-05-22:109   Median : 8.00         2014-06-09:173  
 2014-05-13:106   Mean   :18.41         2014-06-04:154  
 2014-05-21:106   3rd Qu.:10.00         2014-06-06:131  
 2014-05-19:105   Max.   :98.00         (Other)   :139  
 (Other)   :710                         NA's      :359  
          t2.osalusviis t2.evalimiste.usaldus
 E-hääletas      :141   Min.   : 0.00        
 Ei valinud      :243   1st Qu.: 5.00        
 Hääletas paberil:317   Median : 8.00        
 NA's            :660   Mean   :16.48        
                        3rd Qu.:10.00        
                        Max.   :97.00        
                        NA's   :359

Vaatame näiteks, kas kurikuulus „Haldermanni case“ avaldas mõju e-valimiste usaldusele (vt näit Eksperdid: Eesti e-valimised on nii ebaturvalised, et tuleks kohe ära lõpetada). Antud asja kajastus levis laiemalt 13. mail selleks ajaks oli paneeluuringu esimeses laines juba piisavalt palju inimesi intervjueeritud, vaatame kas ja kuidas nende usaldus nägi välja enne ja kuidas pärast.

Esimene intervjueerimislaine kestis ka pärast antud juhtumit, mistõttu tuleks sealt välja filtreerida vaatlused pärast 13. maid. Ühtlasi puhastame andmed keeldunud ja puuduvatest väärtustest (kodeeritud 97 ja 98).

enne.parast = paneel %>%
  mutate(t1.date = as.Date(t1.date),
         usaldus.enne = case_when(
           t1.date < as.Date("2014-05-13") & !(t1.evalimiste.usaldus %in% c(97:98)) ~ t1.evalimiste.usaldus),
         usaldus.parast = case_when(
           !(t2.evalimiste.usaldus %in% c(97:98)) ~ t2.evalimiste.usaldus)
         ) %>%
  select(usaldus.enne, usaldus.parast)

summary(enne.parast[c("usaldus.enne", "usaldus.parast")])
  usaldus.enne    usaldus.parast  
 Min.   : 0.000   Min.   : 0.000  
 1st Qu.: 5.000   1st Qu.: 5.000  
 Median : 8.000   Median : 8.000  
 Mean   : 6.767   Mean   : 6.651  
 3rd Qu.:10.000   3rd Qu.: 9.000  
 Max.   :10.000   Max.   :10.000  
 NA's   :1022     NA's   :468

Kuidas muuta andmeid R-is juhul kui tingimusel. Kõige lihtsam variant on kasutada ifelse(...) käsklust, millesse läheb kolm argumenti:

  1. tingimus (loogiline võrdlus);
  2. väärtus tõese tulemuse korral;
  3. väärtus väära tulemuse korral.

Näiteks võtame vektori, mis sisaldab erinevaid taimi ja kodeerime selle ümber puu- ja köögiviljadeks:

taimed = c("apelsin", "kurk", "õun", "porgand", "peet")
ifelse(taimed %in% c("apelsin", "õun"), "puuvili", "köögivili")
[1] "puuvili"   "köögivili" "puuvili"   "köögivili" "köögivili"

Mida aga teha juhul, kui gruppe on üle kahe? If else puhul tuleks lisada iga järgnev kriteerium väära väärtussesse:

taimed = c("sibul", "till", "küüslauk", "kapsas", "piparmünt")
ifelse(taimed %in% c("sibul", "küüslauk"), "sibulköögivili",
       ifelse(taimed %in% c("till", "piparmünt"), "maitseköögivili", "lehtköögivili"))
[1] "sibulköögivili"  "maitseköögivili" "sibulköögivili"  "lehtköögivili"  
[5] "maitseköögivili"

Gruppide kasvades pole if else enam seega sobilik vahend. Paindlikum on kasutada paketist dplyr tulenevat funktsiooni case_when(), millesse tuleb sisestada vaid kaks argumenti:

  1. tingimus;
  2. väärtus tõese tulemuse korral.

Argumendid tuleb eraldada sümboliga ~ ja erinevad argumendi read komadega. Väärad tulemused märgitakse vaikimise puuduvateks väärtusteks, kuid ka neile võib anda väärtuse, mispuhul oleks tingimuseks lihtsalt TRUE:

case_when(taimed %in% c("sibul", "küüslauk") ~ "sibulköögivili",
          taimed %in% c("till", "piparmünt") ~ "maitseköögivili",
          TRUE ~ "lehtköögivili")
[1] "sibulköögivili"  "maitseköögivili" "sibulköögivili"  "lehtköögivili"  
[5] "maitseköögivili"

case_when argumendi ridu saab esitada lõpmatu arv. Rohkem infot: case_when.

Nüüd saame teha kahe sõltuva valimi t-testi.

t.test(enne.parast$usaldus.enne, enne.parast$usaldus.parast, paired=T)

    Paired t-test

data:  enne.parast$usaldus.enne and enne.parast$usaldus.parast
t = 0.62962, df = 250, p-value = 0.5295
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1865259  0.3618247
sample estimates:
mean of the differences 
              0.0876494

4.2 Anova

Eelmisest ülesandest oli näha, et statistiliselt olulist seost polnud. Kuid et tulemuses päris kindel olla, tuleks võtta arvesse veel teisigi võimalusi:

  1. efekt võis olla lühiajaline, mis tähendab, et selle mõju lahtus järeluuringu ajaks;
  2. efekt võis sõltuda lisatingimustest, mida eelmises näites arvesse ei võetud.

Kontrollime neid variante! Esmalt vaatame järgi, kas efekt oli lühiajaline. Üks viis selleks on kollabeerida andmed intervjueerimise päevade kaupa ja võrrelda muutust ajas.

paneel %>%
  mutate(usaldus = case_when(!(t1.evalimiste.usaldus %in% 97:98) ~ t1.evalimiste.usaldus)) %>%
  group_by(t1.date) %>%
  summarise(keskmine = mean(usaldus, na.rm=T),
            sd = sd(usaldus, na.rm=T)) %>%
  ggplot() +
  geom_line(aes(as.Date(t1.date), keskmine, color="Keskmine")) +
  geom_line(aes(as.Date(t1.date), sd, color="Standardhälve")) +
  scale_x_date(date_breaks="2 days") +
  labs(x="Kuupäev", y="Usaldus (0-10)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=30,vjust=0.6))

Tõepoolest on näha, et usaldus langeb 13. ja 14. mai, kuid seejärel taastub.

Teiseks uurime, kas erinevus võib ilmneda mingite kolmandate tunnuste alusel. Andmestikus on meil selleks tunnus osalusviisi kohta.

summary(paneel$t2.osalusviis)
      E-hääletas       Ei valinud Hääletas paberil             NA's 
             141              243              317              660

Seega tuleks teha juba varasem enne-pärast võrdlus antud gruppide lõikes. Selleks kasutame ANOVA-t, R-is aov(...):

anova.data = paneel %>%
  mutate(t1.date = as.Date(t1.date),
         usaldus.enne = case_when(
           t1.date < as.Date("2014-05-13") & !(t1.evalimiste.usaldus %in% c(97:98)) ~ t1.evalimiste.usaldus),
         usaldus.parast = case_when(
           !(t2.evalimiste.usaldus %in% c(97:98)) ~ t2.evalimiste.usaldus)
         ) %>%
  select(usaldus.enne, usaldus.parast, t2.osalusviis)

aov(usaldus.enne-usaldus.parast ~ t2.osalusviis, data=anova.data) %>% summary()
               Df Sum Sq Mean Sq F value Pr(>F)  
t2.osalusviis   2   35.2  17.624   3.743 0.0257 *
Residuals     170  800.4   4.708                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1188 observations deleted due to missingness

Ka siin selgub, et esineb statistiliselt oluline seos. Paraku ei anna anova otsest aimdust selle kohta, milliste gruppide vahel seos on. Hea on seega vaadata jaotusi graafiliselt, näiteks karpdiagrammi abil.

anova.data %>%
  na.omit() %>%
  ggplot() +
  geom_boxplot(aes(t2.osalusviis, usaldus.enne-usaldus.parast)) +
  theme_bw()

Pilt viitab justkui sellele, et paberhääletajatel jaotusid erinevused võrdselt ümber nulli, samas e-hääletajate usaldus oli langenud ja mittehääletajatel tõusnud.

4.3 Efekti suurus (Coheni d)

Tihti huvitab meid statistikas efekti suurus. St, kui suurt mõju üks tunnus teisele avaldab. Keskmiste võrdluse puhul loetakse efekti suuruseks erinevust erinevust kahe keskmise vahel, jagatuna andmete standardhälvega:

\(d=\frac{\bar{x_1}-\bar{x_2}}{s}\)

Antud efekti suurust nimetatakse Coheni d-ks. Coheni d leidmiseks kasutame paketti effsize ja funktsiooni cohen.d(...), mis jäljendab oma loogikalt t-testi t.test(...). (Antud funktsioon kasutab siiski eelnevalt toodud valemi robustsemat edasiarendust, kuid selle loogika on siiski sama!)

df %>%
  filter(party %in% c("kesk", "ref")) %>%
  cohen.d(support ~ party, data=.)

Cohen's d

d estimate: 0.536684 (medium)
95 percent confidence interval:
inf sup 
NaN NaN

5 Visualiseerimine

5.0.1 Paketid

library(tidyverse)

5.0.2 Andmestik

df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")
paneel = read.csv("https://martenveskimae.github.io/stat-mod/ep2014_paneel.csv")

5.1 Tulpdiagramm

Tulpdiagramm on sobilik juhul, kui soovime võrrelda kategoorilise tunnuse gruppe ühe või väikese arvu tunnuste lõikes. Juhul, kui tunnuseid on mitu, võib tunnuste väärtused asetada üksteise otsa või üksteise kõrvale (lähemalt all). Antud diagrammi eelis seisneb võimaluses hinnata lihtsa vaevaga tulpade erinevusi tulpade kõrgustes (sellest tulenevalt tuleks võrreldavate tunnuste arv hoida pigem madalal).

Püüame näiteks graafida erakondade keskmised toetusnumbrid aastate lõikes. Tulpdiagrammi loomiseks on funktsioon geom_col(...). Teeme tulpdiagrammi, kus x-teljel on aastad, y-teljel iga erakondade keskmine toetus ning tulba värv on seotud erakonnaga.

df %>%
  group_by(party, year) %>%
  summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
  ggplot() +
  geom_col(aes(year, keskmine.toetus, fill=party))

Antud juhul käsitletakse x-teljel olevat tunnust year arvtunnusena, mistõttu on ggplot ise otsustanud, millise vahega väärtusi joonisel kuvada. Praegu võiks kuvada kõik aastad. Kõige lihtsam variant on kasutada käsklust as.factor(year), mis muudab tunnuse kategooriliseks (võib kasutada nii ggploti sees, kui ka enne seda).

df %>%
  group_by(party, year = as.factor(year)) %>%
  summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
  ggplot() +
  geom_col(aes(year, keskmine.toetus, fill=party))

Joonisel on keeruline võrrelda tulpade keskel olevaid erakondi. Parem oleks asetada erinevad erakonnad üksteise kõrvale. Selleks kasutame argumenti position: geom_col(position="dodge").

df %>%
  group_by(party, year = as.factor(year)) %>%
  summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
  ggplot() +
  geom_col(aes(year, keskmine.toetus, fill=party), position="dodge")

Parem! Järgmiseks oleks hea muuta joonise värve üksteisest selgemalt eristuvaks. Muudame ära ka legendi pealkirja. Selleks peaksime muutma fill skaala argumente scale_fill_brewer(palette="Accent", name="Erakonnad") (värvide käsitsi sisestamiseks tuleks kasutada scale_fill_manual(values=c(...))).

R kasutab värvide genereerimiseks Color Brewerit (pakett RColorBrewer). Rohkem infot: R Color Cheatsheet

df %>%
  group_by(party, year = as.factor(year)) %>%
  summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
  ggplot() +
  geom_col(aes(year, keskmine.toetus, fill=party), position="dodge") +
  scale_fill_brewer(palette="Accent", name="Erakonnad")

Veel tuleks muuta y-telje skaala protsentideks, kasutades selleks paketti scales scale_y_continuous(labels=scales::percent). Muudame ära ka telgede nimetused ning lisame pealkirja ja viite käsklusega labs(x="Aasta", y="Toetus",title="Erakondade keskmine toetus",caption="Kantar Emor küsitlusandmetel").

df %>%
  group_by(party, year = as.factor(year)) %>%
  summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
  ggplot() +
  geom_col(aes(year, keskmine.toetus, fill=party), position="dodge") +
  scale_fill_brewer(palette="Accent", name="Erakonnad") +
  scale_y_continuous(labels=scales::percent) +
  labs(x="Aasta", y="Toetus",title="Erakondade keskmine toetus",caption="Kantar Emor küsitlusandmetel")

Viimaks muudame joonise tausta valgeks juba tuttava käsklusega theme_bw().

df %>%
  group_by(party, year = as.factor(year)) %>%
  summarise(keskmine.toetus = mean(support, na.rm=T)) %>%
  ggplot() +
  geom_col(aes(year, keskmine.toetus, fill=party), position="dodge") +
  scale_fill_brewer(palette="Accent", name="Erakonnad") +
  scale_y_continuous(labels=scales::percent) +
  labs(x="Aasta", y="Toetus",title="Erakondade keskmine toetus",caption="Kantar Emor küsitlusandmetel") +
  theme_bw()

5.2 Histogramm

Histogramm on sobilik arvtunnuse jaotuse visualiseerimiseks. Histogrammi kuvamiseks on ggplotis käsklus geom_histogramm(...). Histogrammi y-telg on vaikimisi sagedus, ent selle võib kirjutada soovi korral üle. Histogrammil on võimalik seadistada vahemike laiust või hulka (argumendid bin ja binwidth).

df %>%
  ggplot() +
  geom_histogram(aes(support), binwidth=0.02) +
  scale_x_continuous(labels=scales::percent) +
  labs(x="Toetus", y="Sagedus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
  theme_bw()

Iga erakonna kohta eraldi joonise kuvamiseks võime kasutada käsklusi facet_wrap(~party) või facet_grid(~party).

df %>%
  ggplot() +
  geom_histogram(aes(support), binwidth=0.02) +
  facet_wrap(~party) +
  scale_x_continuous(labels=scales::percent) +
  labs(x="Toetus", y="Sagedus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
  theme_bw()

5.3 Karpdiagramm ja viiuldiagramm

Arvtunnuse jaotuse näitamiseks on võimalik kasutada lisaks histogrammile karpdiagrammi geom_boxplot(...) ja viiuldiagrammi geom_violin(...). Karpdiagrammile võib hajutatult kuvada ka üksikud vaatlused käsklusega geom_jitter(...) (mõistlik on muuta kuvatud vaatlused läbipaistvaks argumendiga alpha ning koondada need kitsamale alale argumendiga width).

df %>%
  ggplot() +
  geom_boxplot(aes(party,support)) +
  geom_jitter(aes(party,support), alpha=0.2, width=0.2) +
  scale_y_continuous(labels=scales::percent) +
  labs(x="Erakond", y="Toetus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
  theme_bw()

Viiuldiagramm annab veelgi täpsema ülevaate jaotusest. Viiuldiagrammile võib lisada kvantiilide jooned argumendiga draw_quantiles.

df %>%
  ggplot() +
  geom_violin(aes(party,support), draw_quantiles=c(.25,.5,.75)) +
  scale_y_continuous(labels=scales::percent) +
  labs(x="Erakond", y="Toetus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
  theme_bw()

5.4 Hajuvusdiagramm

Hajuvusdiagrammi kasutatakse kahe arvtunnuse omavahelise seose visualiseerimiseks.

df %>%
  ggplot() +
  geom_point(aes(support,gdpgr_yearly)) +
  scale_x_continuous(labels=scales::percent) +
  facet_wrap(~party) +
  labs(x="Toetus", y="SKT aastane muutus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
  theme_bw()

5.5 Joondiagramm

Joondiagramm, nagu hajuvusdiagrammgi, on sobilik kahe arvtunnuse omavahelise seose visualiseerimiseks. Joondiagrammi kasutatakse rohkem siis, kui üks tunnustest on seotud ajaga.

df %>%
  ggplot() +
  geom_line(aes(as.Date(date),support)) +
  scale_y_continuous(labels=scales::percent) +
  facet_wrap(~party) +
  labs(x="", y="Toetus", title="Erakondade toetus",caption="Kantar Emor küsitlusandmetel") +
  theme_bw()

5.5 Likerti skaala

Küsitlustes kasutatakse tihti likerti skaalasid, mida visualiseeritakse enamasti tulpdiagrammi abil (skaala on protsentides).

Kasutame näiteks EP 2014 andmestikust usalduse tunnust.

library(reshape2)
usaldus = paneel %>%
  mutate(enne = case_when(!(t1.evalimiste.usaldus %in% 97:98) ~ t1.evalimiste.usaldus),
         parast = case_when(!(t2.evalimiste.usaldus %in% 97:98) ~ t2.evalimiste.usaldus)) %>%
  select(enne,parast) %>%
  melt() %>%
  na.omit() %>%
  group_by(variable,value = as.factor(value)) %>%
  summarise(count = n()) %>%
  mutate(pc = count/sum(count[variable==variable]))

usaldus
# A tibble: 22 x 4
# Groups:   variable [2]
   variable  value count         pc
     <fctr> <fctr> <int>      <dbl>
 1     enne      0   185 0.15651438
 2     enne      1    17 0.01438240
 3     enne      2    29 0.02453469
 4     enne      3    44 0.03722504
 5     enne      4    32 0.02707276
 6     enne      5    85 0.07191201
 7     enne      6    42 0.03553299
 8     enne      7    88 0.07445008
 9     enne      8   185 0.15651438
10     enne      9   176 0.14890017
# ... with 12 more rows

Graafime tulemused tulpdiagrammina.

usaldus %>%
  ggplot() +
  geom_col(aes(variable,pc,fill=value))

Antud joonist saab muuta loetavamaks vahetades värviskaalat, pöörates joonise telgi, muutes telgede nimesid jne.

usaldus %>%
  ggplot() +
  geom_col(aes(variable,pc,fill=reorder(value,desc(value))),width=0.7) +
  scale_fill_brewer(name="Usaldus",palette="Paired") +
  scale_y_continuous(labels=scales::percent) +
  scale_x_discrete(labels=c("Enne","Pärast")) +
  labs(x="", y="Osakaal", title="E-valimiste usaldus",caption="EP 2014 valimiste küsitlusandmed") +
  coord_flip() +
  theme_minimal()

Alternatiiv on kasutada spetsiaalset likerti skaalade visualiseerimiseks mõeldud funktsiooni likert(...) paketist HH.

usaldus.2 = usaldus %>%
  dcast(variable ~ value, value.var="count")
rownames(usaldus.2) = c("Enne", "Pärast")

HH::likert(usaldus.2,
           as.percent=T,
           ylab.right="Vastanute arv",
           xlab="Osakaal",
           main="E-valimiste usaldus",
           sub="")

6 Korrelatsioon

6.0.1 Paketid

library(tidyverse)
library(psych)
library(PerformanceAnalytics)

6.0.2 Andmestik

df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")
paneel = read.csv("https://martenveskimae.github.io/stat-mod/ep2014_paneel.csv")

6.1 Pearsoni korrelatsioonikordaja

Vaatame lähemalt kahe tunnuse vahelist korrelatsiooni, antud juhul Pearsoni korrelatsiooni. Pearsoni korrelatsioon sobib juhul, kui mõlema tunnuse puhul on tegu arvtunnustega.

Pearsoni korrelatsioon on iseloomult lineaarne, st kahe tunnuse vahele üritatakse sobitada lineaarset seost. Matemaatiliselt on korrelatsioon kahe tunnuse kovariatsioon, skaleerituna nende standardhälbega:

\(r=\frac{\sum\limits_{i=1}^n(x_{i}-\bar{x})(y_{i}-\bar{y})}{(n-1)s_{x}s_{y}}\)

Et anda hinnang korrelatsioonikordaja statistilisele olulisusele, võime arvutada korrelatsiooni t-statistiku väärtuse ja võrrelda seda teoreetilise t-jaotusega.

\(t=\frac{r\sqrt{n-2}}{\sqrt{1-r^2}}\)

Arvutame käsitsi kahe juhusliku vektori korrelatsiooni (eeldatavalt on korrelatsioon juhuslikkuse tõttu 0). Esmalt tekitame kaks vektorit ning vaatame võimalikku visuaalselt.

n = 100
x = rnorm(n)
y = rnorm(n)
plot(x,y)

Kui kahe vektori vahel on tõepoolest seos, siis on see juhuslik. Hindame seost korrelatsiooni abil.

xbar = mean(x)
ybar = mean(y)
sx = sd(x)
sy = sd(y)

r = sum((x-xbar)*(y-ybar)) / ((n-1)*sx*sy)
r
[1] 0.1136141

Kas tulemus on statistiliselt oluline?

t = (r*sqrt(n-2)) / sqrt(1-r^2) # t-väärtus
p = 2*pt(-abs(t),df=199) # p-väärtus

# Olulisuse visualiseerimine (sama, mis t-testi puhul)
tdistx = seq(-3,3,length=300)
tdisty = dt(tdistx,df=199)
a95 = qt(c(.025, .975), df=199)
ggplot() +
  geom_line(aes(tdistx,tdisty)) +
  geom_vline(xintercept = a95[1]) +
  geom_vline(xintercept = a95[2]) +
  geom_vline(xintercept = t, color="red") +
  theme_bw()

Arvutuste tulemusel on:

  • korrelatsioon: 0.1136141
  • t-skoor: 1.1320519
  • p-väärtus: 0.2589748

Õnneks teeb kogu töö ära funktsioon corr.test(...) paketist psych:

corr.test(cbind(x,y))
Call:corr.test(x = cbind(x, y))
Correlation matrix 
     x    y
x 1.00 0.11
y 0.11 1.00
Sample Size 
[1] 100
Probability values (Entries above the diagonal are adjusted for multiple tests.) 
     x    y
x 0.00 0.26
y 0.26 0.00

 To see confidence intervals of the correlations, print with the short=FALSE option

Tulemus on sama!

Majandusliku hääletamise teooriast teame, et hääletamine on tihti seotud majanduse käekäiguga - kui majandusel läheb hästi, on koalitsiooni toetus kõrge, kui halvasti, on opositsiooni toetus kõrge. Vaatame näiteks, kuidas on seotud Reformierakonna toetus tööpuuduse, SKT kasvu ja inflatsiooniga.

ref = df %>%
  filter(party == "ref") %>%
  select(support, unemp, gdpgr_yearly, inflats_m) %>%
  na.omit()

cor(ref)
                support      unemp gdpgr_yearly  inflats_m
support       1.0000000 -0.0712523    0.3312793  0.3958776
unemp        -0.0712523  1.0000000   -0.1915624 -0.2088902
gdpgr_yearly  0.3312793 -0.1915624    1.0000000  0.3311977
inflats_m     0.3958776 -0.2088902    0.3311977  1.0000000

Korrelatsioonide statistilise olulisuse jaoks tuleks kasutada corr.test(...) funktsiooni paketis psych.

corr.test(ref)
Call:corr.test(x = ref)
Correlation matrix 
             support unemp gdpgr_yearly inflats_m
support         1.00 -0.07         0.33      0.40
unemp          -0.07  1.00        -0.19     -0.21
gdpgr_yearly    0.33 -0.19         1.00      0.33
inflats_m       0.40 -0.21         0.33      1.00
Sample Size 
[1] 93
Probability values (Entries above the diagonal are adjusted for multiple tests.) 
             support unemp gdpgr_yearly inflats_m
support          0.0  0.50         0.01      0.00
unemp            0.5  0.00         0.13      0.13
gdpgr_yearly     0.0  0.07         0.00      0.01
inflats_m        0.0  0.04         0.00      0.00

 To see confidence intervals of the correlations, print with the short=FALSE option

Vaikimisi annab antud funktsioon nii korrelatsioonimaatriksi kui p-väärtused. Kui soovida vaid ühte, tuleks käskluse lõppu lisada kas $r või $p korrelatsioonide või p-väärtuste saamiseks.

corr.test(ref)$r
                support      unemp gdpgr_yearly  inflats_m
support       1.0000000 -0.0712523    0.3312793  0.3958776
unemp        -0.0712523  1.0000000   -0.1915624 -0.2088902
gdpgr_yearly  0.3312793 -0.1915624    1.0000000  0.3311977
inflats_m     0.3958776 -0.2088902    0.3311977  1.0000000

Alati on mõistlik hinnata pilti graafiliselt. Korrelatsiooni puhul kasutatakse tihti spetsiaalseid korrelatsiooni visualiseeringuid, mida on keeruline ggplotis luua. Küll aga on mitmeid pakette, mis antud probleemi lahendavad. Näiteks võib kasutada funktsiooni chart.Correlation(...) paketist PerformanceAnalytics.

chart.Correlation(ref)

6.2 Spearmani korrelatsioon

Juhul, kui tunnused on järjestikused, on mõistlik kasutada Spearmani korrelatsiooni, mis ei kasuta korrelatsiooni arvutamiseks algtunnused, vaid annab neile astmeväärtuse (rank). Spearmani korrelatsiooni läbiviimiseks võib kasutada juba tuttavat EP 2014 valimiste andmestikku ning usalduse tunnuseid.

usaldus = paneel %>%
  mutate(enne = case_when(!(t1.evalimiste.usaldus %in% 97:98) ~ t1.evalimiste.usaldus),
         parast = case_when(!(t2.evalimiste.usaldus %in% 97:98) ~ t2.evalimiste.usaldus)) %>%
  select(enne,parast) %>%
  na.omit()

corr.test(usaldus, method="spearman")
Call:corr.test(x = usaldus, method = "spearman")
Correlation matrix 
       enne parast
enne   1.00   0.75
parast 0.75   1.00
Sample Size 
[1] 835
Probability values (Entries above the diagonal are adjusted for multiple tests.) 
       enne parast
enne      0      0
parast    0      0

 To see confidence intervals of the correlations, print with the short=FALSE option
chart.Correlation(usaldus, method="spearman")

7 Lineaarne regressioon

7.0.1 Paketid

library(tidyverse)
library(reshape2)

7.0.2 Andmestik

df = read.csv("https://martenveskimae.github.io/stat-mod/party_dataset.csv")

7.1 Mudeli spetsifitseerimine R-is

Mitmed mudelid R-is, sh lineaarse regressiooni mudel, küsivad sisendit valemi kujul. Meenutame näiteks lineaarse regressiooni valemit:

\(\hat{y} = \alpha + \beta X + \epsilon\)

kus:

  • \(\hat{y}\) on ennustatav sõltuva tunnuse väärtus;
  • \(\alpha\) on sõltuva tunnuse väärtus, kui sõltumatud tunnused on null. Tihti tähistatud ka kui \(\beta_0\);
  • \(\beta\) on sõltumatute tunnuste hinnanguline kordaja, teisisõnu hinnanguline mõju suurus;
  • \(X\) on sõltumatud tunnused;
  • \(\epsilon\) on viga, mida enamasti käsitletakse juhusliku “mürana” andmetes.

Lineaarse regressiooni jooksutamiseks oleks vaja sisendiks anda nii y kui ka x tunnused, misjärel arvutatakse välja \(\alpha\) ja \(\beta\) väärtused.

7.2 Erakondade toetus ja makromajanduslikud näitajad

Jätkame siinkohal korrelatsioonanalüüsis uuritud tunnustega - erakonna toetus ja töötus. Näitena kasutame taas Reformierakonda. Valemi kujul näeks meie mudel välja järgnev:

\(toetus_{ref} = \beta_0 + \beta_1 töötus + \epsilon\)

Mudeli spetsifitseerimisel näeks see välja toetus ~ töötus, kus ~ sümbol eristab sõltuvat ja sõltumatuid tunnuseid. Lineaarse regressiooni jooksutamiseks kasutame funktsiooni lm(...), kuhu sisendiks läheb valem ning andmestiku nimetus (data=...).

Filtreerime andmestikust Reformierakonnaga seotud vaatlused ning kasutame support ja unemp tunnuseid. Enne analüüsini jõudmist viime veel support tunnuse 0 ja 100 vahele, et ühtlustada erinevate tunnuste skaalad ning lihtsustada hilisemat tulemuste tõlgendamist. Lisaks arvutame aastase inflatsiooni, mida kasutame mitmese lineaarse regressiooni näites.

df = df %>%
  mutate(support=support*100, 
         inflats_a =((inflats_m_basevalue / inflats_prevyear_base) - 1)*100)

ref = df %>% filter(party=="ref")

ref_mudel = lm(support ~ unemp, data=ref)
ref_mudel

Call:
lm(formula = support ~ unemp, data = ref)

Coefficients:
(Intercept)        unemp  
    27.3345       0.2048

Mudeli jooksutamisel kuvatakse \(\beta_0\) (Intercept) ja \(\beta_1\) (unemp), antud juhul väärtustega 27.335 ja 0.205. Kuna tunnused on nüüd samal protsentuaalsel skaalal, siis saame järeldada, et töötuse 1% suurenedes kasvab Reformierakonna toetus 0.2%. Arvestades, et Reformierakond on olnud pikalt valitsuses, siis võiks eeldada hoopis negatiivset seost - töötuse kasvades langeb valitseva erakonna toetus. Mudeliga tuleks seega veel tööd teha, ent esmalt uurime teisi mudeliga seotud statistikuid.

summary(ref_mudel)

Call:
lm(formula = support ~ unemp, data = ref)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.6659  -5.0756   0.2317   3.7196  16.6618 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  27.3345     1.6422  16.645   <2e-16 ***
unemp         0.2048     0.1615   1.268    0.208    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.772 on 107 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared:  0.01481,   Adjusted R-squared:  0.005599 
F-statistic: 1.608 on 1 and 107 DF,  p-value: 0.2075

summary(...) funktsioon annab meile lisaks \(\beta\) väärtuste standard vead ja p-väärtused, ning \(R^2\) väärtused. Näeme, et töötuse tunnus ei ole statistiliselt oluline (\(\beta_1\) ei erine oluliselt nullist) ning mudeli seletusvõime on vaid 1.5% (Multiple R-squared). Ühtlasi näeme viimaselt realt, et kõiki sõltumatuid tunnuseid arvesse võttes pole mudel statistiliselt oluline (p-väärtus: 0.208). Kuna antud juhul kaasasime vaid ühe sõltumatu tunnuse, ei paku see uut teadmist, ent keerulisemate mudelite puhul, kus on mitmeid sõltumatuid tunnuseid, tasub seda jälgida.

Regressioonanalüüsis hinnatakse \(\beta\) väärtuste statistilist olulisust t-testi abil, hinnates, kas väärtus on oluliselt erinev nullist. Sellest tulenevalt, mida lähemal on \(\beta\) nullile, seda suurem kipub olema tema p-väärtus. Seega, väikeste \(\beta\) väärtuste puhul tuleks pidada meeles, et statistiliselt olulise tulemuse saamise tõenäosus on madalam.

Alljärgnevalt viime läbi sumulatsiooni ning uurime \(\beta_0\) ja \(\beta_1\) varieeruvust erinevate valimi suuruste puhul. Hindame seejuures tüüp I viga ning näeme, et see vastab valitud usaldusnivoole.

Populatsiooni parameetrite suurusi võib muuta, et piltlikustada eelpool kirjeldatud nähtust.

N = 100000          # Populatsiooni suurus
n = 1000            # Valimi ülemine piir
x = rnorm(N,50,15)  # X
e = rnorm(N,0,10)   # Juhuslik viga
b0 = 2.5            # beta 0
b1 = 5              # beta 1
y = b0 + (b1*x) + e # Y

# Loome andmetabeli x ja y väärtustega
pop = data.frame(x = x, y = y)

# Loome funktsiooni, mis väljastab regressiooni tulemused vastavalt valimi suurusele
reg = function(data,n){
  df = data[sample(1:nrow(data), n),]
  m = lm(y~x,data=df)
  data.frame(b0 = m$coefficients[[1]],
             b0l = confint(m)[1,1],
             b0u = confint(m)[1,2],
             b0p = summary(m)$coefficients[,4][[1]],
             b0se = summary(m)$coefficients[,2][[1]],
             b1 = m$coefficients[[2]],
             b1l = confint(m)[2,1],
             b1u = confint(m)[2,2],
             b1p = summary(m)$coefficients[,4][[2]],
             b1se = summary(m)$coefficients[,2][[2]],
             n=n)
}

# Jooksutame regressiooni erinevate valimi suurustega
r = lapply(10:n,function(x)reg(pop,x)) %>% do.call(rbind.data.frame,.)

# Hindame, kas saadud tulemustes esines tüüp I viga (false positive)
b0_pop = reg(pop,N)$b0
b1_pop = reg(pop,N)$b1

r = r %>%
  mutate(b0_typeI = ifelse((b0l > b0_pop | b0u < b0_pop) & b0p <= .05, TRUE, FALSE),
         b1_typeI = ifelse((b1l > b1_pop | b1u < b1_pop) & b1p <= .05, TRUE, FALSE))

# Graafime tulemused
r %>%
  select(n,b0_typeI,b1_typeI,b0,b1) %>%
  melt(c("n","b0_typeI","b1_typeI")) %>%
  melt(c("n","variable","value"),variable.name="TF",value.name="type_I") %>%
  ggplot() +
  geom_point(aes(n,value,color=type_I),shape=1,alpha=.5) +
  geom_hline(data=data.frame(variable=c("b0","b1"),
                             value=c(b0_pop,b1_pop)),
             aes(yintercept=value)) +
  scale_color_manual(values=c("steelblue","red")) +
  facet_wrap(~variable) +
  theme_bw()

r %>%
  ggplot() +
  geom_point(aes(n,b0p,color=b0_typeI),shape=1) +
  geom_hline(yintercept=0.05) +
  scale_color_manual(values=c("steelblue","red")) +
  theme_bw()

# Arvutame, mitmel protsendil juhtudest esines tüüp I viga (eeldada võiks ligemale 5%)
sum(r$b0_typeI==TRUE)/nrow(r)
[1] 0.02825429
sum(r$b1_typeI==TRUE)/nrow(r)
[1] 0.06559031

Tihti kasutatakse mudeli hindamiseks veel ruutkeskmist hälvet (mean squared error ehk MSE), mis väljendab, kui hästi suudab mudel ennustada sõltuva tunnuse väärtuseid. MSE leidmiseks tuleb lahutada õigetest väärtustest ennustatud väärtused (ehk arvutada ennustusjäägid), võtta need ruutu ning leida saadud jada keskmine. Seda on lihtne teha, kuivõrd mudel pakub ennustusjääke $residuals abil.

mean(ref_mudel$residuals^2)
[1] 45.01441

Kui soovime näha mudeli ennustatud väärtusi, siis on selleks funktsioon predict(...), kuhu läheb argumendiks mudel ning soovi korral andmestik uute andmetega.

yhat = predict(ref_mudel)
head(yhat)
       1        2        3        4        5        6 
28.42012 28.33819 28.33819 28.33819 28.17432 28.17432

Sarnase mudeli võib teha ka teiste erakondade kohta. Võtame näiteks Keskerakonna.

kesk = df %>% filter(party=="kesk")

kesk_mudel = lm(support ~ unemp, data=kesk)
summary(kesk_mudel)

Call:
lm(formula = support ~ unemp, data = kesk)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.1163 -2.3599 -0.3058  2.1375  9.5876 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 24.53591    0.82636  29.692   <2e-16 ***
unemp        0.11844    0.08128   1.457    0.148    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.408 on 107 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared:  0.01946,   Adjusted R-squared:  0.0103 
F-statistic: 2.124 on 1 and 107 DF,  p-value: 0.148

Tulemus on suuresti sama.

Mudeli tulemuste mugavaks esitamiseks on loodud pakett stargazer. Stargazer võimaldab esitada nii lineaarse regressiooni kui mitmete teiste mudelite tulemusi kujul, mis on aktsepteeritav enamikes teadusajakirjades. Stargazeri võimalustega tutvumiseks saab infot järgnevalt lingilt: Stargazer cheatsheet.

library(stargazer)

# type võimalikud väärtused on "text", "html" ja "latex", vastavalt soovitud väljundi formaadile (text, html ja pdf)
stargazer(ref_mudel, kesk_mudel,
          title="Linear regression table",
          covariate.labels="Unemployment",
          dep.var.labels="Support",
          column.labels=c("Ref", "Kesk"),
          model.numbers=F,
          no.space=T,
          star.cutoffs=c(0.05,0.01,0.001),
          type="html")
Linear regression table
Dependent variable:
Support
Ref Kesk
Unemployment 0.205 0.118
(0.162) (0.081)
Constant 27.335*** 24.536***
(1.642) (0.826)
Observations 109 109
R2 0.015 0.019
Adjusted R2 0.006 0.010
Residual Std. Error (df = 107) 6.772 3.408
F Statistic (df = 1; 107) 1.608 2.124
Note: p<0.05; p<0.01; p<0.001

7.3 Mitmene lineaarne regressioon

Lineaarset regressiooni kasutades ei pea piirduma vaid ühe tunnuse lisamisega, vaid see võimaldab hinnata mitmeid tunnuseid korraga (oluline eelis t-testi ja korrelatsiooni ees). Lisame varasemale mudelile SKT kasvu ja inflatsiooni tunnused:

\(toetus_{ref} = \beta_0 + \beta_1 töötus + \beta_2 SKT + \beta_3 inflatsioon + \epsilon\)

mis näeks R-i valemina välja support ~ unemp + gdpgr_yearly + inflats_m. Lihtsustamise mõttes piirame vaatlused eelviimase valitsustsükliga, et kontrollida koalitsiooni/opositsiooni mõju ning jooksutame uue mudeli lisatunnustega.

ref = ref %>% filter(year>=2011, year<2015)
  
mlr_mudel = lm(support ~ unemp + gdpgr_yearly + inflats_a, data=ref)
summary(mlr_mudel)

Call:
lm(formula = support ~ unemp + gdpgr_yearly + inflats_a, data = ref)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2222 -2.3461  0.6219  1.5866 11.2140 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   29.9480     3.8262   7.827 8.35e-10 ***
unemp         -0.8965     0.5845  -1.534    0.132    
gdpgr_yearly   1.8274     0.3438   5.315 3.59e-06 ***
inflats_a     -0.3587     0.4936  -0.727    0.471    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.535 on 43 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4927,    Adjusted R-squared:  0.4573 
F-statistic: 13.92 on 3 and 43 DF,  p-value: 1.761e-06

Tulemused näitavad, et SKT kasv on statistiliselt oluline ning märkimisväärse mõjuga. Samuti kasvas mudeli seletusvõime 46%ni (Adjusted R-squared) ning mudel ise on statistiliselt oluline.

mean(mlr_mudel$residuals^2)
[1] 11.43551

Mudeli ruutkeskmine hälve on samuti oluliselt vähenenud võrreldes eelnevaga.

Kui võtta eesmärgiks saada võimalikult hea mudeli sobivus, võib kaaluda sõltumatute tunnuste transformeerimist, eeldusel, et see on teoreetiliselt põhjendatud. Seega tuleks küsida, milline funktsionaalne vorm võiks antud sõltumatutel tunnustel peale lineaarse vormi olla. Visualiseerime seosed.

ref %>%
  melt("support",c("unemp","gdpgr_yearly","inflats_a")) %>%
  ggplot() +
  geom_point(aes(value,support)) +
  coord_fixed() +
  facet_wrap(~variable,scales="free",ncol=2)

Antud juhul võib jätta töötuse lineaarseks, SKT kasvu puhul katsetada logaritmimist ning inflatsiooni puhul kasutada ruutseost. Loome uued transformeeritud tunnused, ning vaatame, kas seos erakonna toetusega on muutunud lineaarsemaks.

ref = ref %>%
  mutate(gdpgr_yearly_log = log(gdpgr_yearly),
         inflats_a_squared = inflats_a^2)

ref %>%
  melt("support",c("unemp","gdpgr_yearly_log","inflats_a_squared")) %>%
  ggplot() +
  geom_point(aes(value,support)) +
  coord_fixed() +
  facet_wrap(~variable,scales="free",ncol=2)

Tundub, et on olnud abi ning võime teha uue mudeli transformeeritud tunnustega.

mlr_mudel_tr = lm(support ~ unemp + gdpgr_yearly_log + inflats_a + inflats_a_squared, data=ref)
summary(mlr_mudel_tr)

Call:
lm(formula = support ~ unemp + gdpgr_yearly_log + inflats_a + 
    inflats_a_squared, data = ref)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.9753 -1.9009 -0.3334  1.8826 10.5300 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        29.0877     3.8250   7.605 2.01e-09 ***
unemp              -0.3867     0.5140  -0.752  0.45602    
gdpgr_yearly_log    3.0627     0.9127   3.356  0.00169 ** 
inflats_a          -2.5772     1.2790  -2.015  0.05033 .  
inflats_a_squared   0.5346     0.2694   1.984  0.05379 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.481 on 42 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.5195,    Adjusted R-squared:  0.4738 
F-statistic: 11.35 on 4 and 42 DF,  p-value: 2.461e-06
mean(mlr_mudel_tr$residuals^2)
[1] 10.83034

Mudeli parameetrites on näha väikest arengut. Edasi tuleks anda põhjalikum hinnang mudeli toimimisele.

7.4 Mudeli hindamine

Kõige kiirema ülevaate võimalikest probleemidest saab ennustusjääkide visualiseerimisel.

plot(mlr_mudel_tr, which=1) # Rohkem joonised ennustusjääkidest saab jättes 'which' täpsustamata

Ennustusjääkide puhul soovime näha, et need oleksid jaotunud ühtlase hajususega ümber nulli. Antud juhul on näha, et ennustusjäägid on siiski teatud määral süstemaatiliselt kaldu, mis viitab probleemidele mudeli funktsionaalses vormis. Samas heteroskedastiivsust ei paista eriti olevat.

Mudeli funktsionaalse vormi hindamiseks kasutatakse tihti Ramsey RESET-testi, mis lisab mudelile regressorite ruutseosed ning hindab, kas see on algmudelist parem. Kasutame selleks funktsiooni resettest(...) paketist lmtest.

library(lmtest)
resettest(mlr_mudel, power=2, type="fitted")

    RESET test

data:  mlr_mudel
RESET = 6.2981, df1 = 1, df2 = 42, p-value = 0.01603
resettest(mlr_mudel_tr, power=2, type="fitted")

    RESET test

data:  mlr_mudel_tr
RESET = 2.4592, df1 = 1, df2 = 41, p-value = 0.1245

Antud juhul viitavad tulemused, et ilma transformeeritud tunnusteta mudelil oli halvem funktsionaalne vorm kui transformeeritud tunnustega mudelil. Transformeeritud tunnustega mudel ei muutu oluliselt paremaks, kui lisada ruutseostega tunnused. Samas võib RESET testi tulemus olla oluline kuupseoste (argument power=2:3) korral, kuivõrd visuaalselt oli näha ennustusjääkide teatud süstemaatilist kõikumist ümber nulli.

Heteroskedastiivsuse hindamiseks on populaarsemad White test ja Breusch-Pagan test, millest võib antud juhul katsetada viimasega. Breusch-Pagani testi puhul sobitatakse uus mudel algmudeli ennustusjääkidele ning hinnatakse, kas sealt ilmnevad olulised seosed. lmtest paketist on selleks funktsioon bptest(...).

bptest(mlr_mudel_tr)

    studentized Breusch-Pagan test

data:  mlr_mudel_tr
BP = 5.2184, df = 4, p-value = 0.2656

Kuna antud juhul osutus test statistiliselt ebaoluliseks, võib järeldada, et heteroskedastiivsusega olulisi probleeme pole.

Viimaks võime kontrollida, ega uute tunnuste lisamisel pole mudelis ilmnenud multikollineaarsust. Multikollineaarsuse korral on sõltumatud tunnused on omavahel tugevas korrelatsioonis, mistõttu on raske eristada üksikute tunnuste mõju (korrelatsioonist tingituna, kui üks väärtus kasvab, kasvavad ka teised väärtused ning üksikud mõjud pole enam eristatavad). Sellest tulenevalt pole hinnatud parameetrid (\(\beta\) väärtused) enam usaldatavad. Kuigi teatud korrelatsioon esineb tunnuste vahel peaaegu alati, on oluline, et sõltumatute tunnuste omavaheline korrelatsioon ei ületaks sõltumatu ja sõltuva tunnuse korrelatsiooni.

Regressioonanalüüsi korral on multikollineaarsuse hindamiseks kasutusel VIF (Variance Inflation Factor) test. VIF test hindab, kui suur oli sõltumatu tunnuse \(\beta_k\) variatsioon teiste sõltumatute tunnuste kontekstis (\(Var(\beta_k)\)) ja kui suur isoleeritult (\(Var_{min}(\beta_k)\)). Seejärel jagab ta esimese teisega \(\frac{Var(\beta_k)}{Var_{min}(\beta_k)} = \frac{1}{1-R_k^2}\). Üldjuhul leitakse, et kui tulemus on üle kümne, võib mudelis esineda multikollineaarsuse probleem.

VIF testi funktsiooni vif(...) saame paketist HH. Kuna HH ja tidyiverse pakettides on osa funktsioone sama nimega, siis ei hakka praegu HH-d eraldi aktiveerima, vaid rakendame sealt ainult hetkel soovitud vif funktsiooni.

HH::vif(mlr_mudel_tr)
            unemp  gdpgr_yearly_log         inflats_a inflats_a_squared 
         4.155656          2.647484         20.732871         27.740997

Antud mudeli puhul on probeelm inflatsiooni ja inflatsioon ruudus tunnustega, mis on samas ootuspärane tulemus (võimegi eeldada tugevat seost antud kahe tunnuse vahel)

HH::vif(mlr_mudel)
       unemp gdpgr_yearly    inflats_a 
    5.210994     3.358901     2.994678

Transformeerimata tunnustega mudelis sellist probleemi ei esinenud ning ilmselt oleksid selle tulemused kokkuvõttes usaldusväärsemad.

8 Logistiline regressioon

8.0.1 Paketid

library(tidyverse)
library(stargazer)
library(ResourceSelection)
library(caret)

8.0.2 Andmestik

RK_2015 = read.csv("https://martenveskimae.github.io/stat-mod/RK_2015.csv")

8.1 Üldistatud lineaarne mudel (generalized linear model)

Eelnevalt käsitletud lineaarse regressiooni mudeli üheks peamiseks eelduseks oli, et ennustusjäägid peavad olema normaaljaotuslikud. Tuleb aga välja, et taolist lineaarset mudelit on võimalik üldistada ning eeldada seejuures väga erinevate jaotustega ennustusjääke (ennustusjäägid peavad olema eksponentide perekonnast). See võimaldab ühe tüüpmudeli abil modelleerida näiteks pidevtunnustega andmeid (normaaljaotus), loendusandmeid (Poissooni jaotus), eksponentseoseid, kui ka binaarseid andmeid (Bernoulli jaotus). Mudeli funktsionaalne vorm on järgnev:

\(E(y)=g^{-1}(X\beta)\)

kus:

  • \(E(y)\) on keskmine sõltuvtunnuse väärtus sõltuvalt \(X\) väärtustest;
  • \(X\beta\) X on sõltumatud tunnused korrutatuna hinnatava parameetriga \(\beta\);
  • \(g\) on tõenäosuse funktsioon, mille abil seotakse keskmine sõltuvtunnuse väärtus sõltumatute tunnustega \(X\beta\) (link function)

Antud juhul keskendume jah/ei jaotusele ehk Bernoulli jaotusele, kus

\(g(\mu) = X\beta = ln\big(\frac{\mu}{1-\mu}\big)\) ja

\(\mu = \frac{1}{1+exp(-X\beta)}\)

Logistiline regressioon on olemuselt lineaarne klassifikaator, st ta sobitab kahe klassi vahele lineaarse eralduspiiri. Selle visualiseerimiseks võib teha simulatsiooni, kus võrdleme logistilist regressiooni mõne tuntuma mitte-lineaarse klassifikaatoriga (juhumets ehk random forest ja tugivektormasin ehk support vector machine).

library(randomForest)
library(e1071)

cos_func = function(x,y){
  set.seed(200212)
  x > y + cos(x)*rnorm(1,0,2)
}
pr = function(x) predict(x,test) %>% as.numeric() -1

data = expand.grid(x=seq(0,10,.25),
                   y=seq(0,10,.25)) %>%
  mutate(label = cos_func(x,y))

log_fit = glm(label~.,data=data,family=binomial,control=glm.control(maxit=20))
rf_fit = randomForest(label~.,data=data,mtry=2)
svm_fit = svm(factor(label)~.,data=data,kernel="radial",probability=T,cost=5)
  
data %>%
  mutate(logit_preds = predict(log_fit,data,"response")>.5,
         rf_preds = predict(rf_fit,data,"response")>.5,
         svm_preds = predict(svm_fit,data)) %>%
  reshape2::melt(c("x","y")) %>%
  ggplot() +
  geom_point(aes(x,y,color=value)) +
  coord_fixed() +
  facet_wrap(~variable) +
  theme_bw()

8.2 Valimisosaluse analüüsimine

Valimisosalust mõõdetakse indiviiditasandil enamasti binaarsel kujul - kas vastaja hääletas või mitte. Kui soovime seda modelleerida, siis saame arvutada küsitlusele vastaja tõenäosuse (0-1) kuuluda kas ühte või teise rühma.

Valimisosaluse analüüsiks võtame ette 2015. aasta Riigikogu valimiste küsitluse ning püüame ennustada inimeste valimisosalust sotsiaal-demograafiliste tunnuste abil. Järjestame esmalt ordinaalsed tunnused ja tutvume andmestikuga.

palgad = c("Alla 200 €","200–275 €","276–350 €","351–425 €","426–500 €","501–575 €",
           "576–700 €","701–1000 €","1001–1300 €","1301–1600 €","1601–1900 €",
           "1901–2200 €","2201–2500 €","Üle 2500 €")
haridus = c("Põhiharidus","Keskharidus","Kõrgharidus")
RK_2015$wage = factor(RK_2015$wage,levels=palgad)
RK_2015$education = factor(RK_2015$education,levels=haridus)

summary(RK_2015)
            vote          age          gender          education  
 Ei hääletanud:150   Min.   :18.00   Mees :350   Põhiharidus:118  
 Hääletas     :646   1st Qu.:37.00   Naine:446   Keskharidus:417  
                     Median :53.00               Kõrgharidus:261  
                     Mean   :51.97                                
                     3rd Qu.:67.00                                
                     Max.   :90.00                                
                                                                  
          ethnicity            wage    
 Eestlane      :653   701–1000 € :143  
 Mitte-eestlane:143   576–700 €  :123  
                      276–350 €  : 96  
                      1301–1600 €: 80  
                      1001–1300 €: 78  
                      351–425 €  : 72  
                      (Other)    :204

Sõltuvtunnus on antud hetkel vote, mille puhul näeme, et 646 küsitletutest hääletas ja 150 ei hääletanud. Viime läbi logistilise regressiooni kõikide andmestikus kaasatud tunnustega. Logistilise regressiooni jaoks kasutame funktsiooni glm(...), kus tuleb lisaks valemile ja andmestikule märkida ennustusjääkide perekond family=binominal.

f = formula(vote ~ age + I(age^2) + gender + education + ethnicity + as.numeric(wage))
osalus_mudel = glm(f, data= RK_2015, family=binomial)
summary(osalus_mudel)

Call:
glm(formula = f, family = binomial, data = RK_2015)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3317   0.3497   0.5297   0.6659   1.4727  

Coefficients:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -1.537e+00  7.877e-01  -1.951 0.051055 .  
age                      3.143e-02  3.188e-02   0.986 0.324268    
I(age^2)                -2.829e-05  3.184e-04  -0.089 0.929201    
genderNaine             -8.873e-02  1.950e-01  -0.455 0.649071    
educationKeskharidus     8.422e-01  2.441e-01   3.449 0.000562 ***
educationKõrgharidus     1.626e+00  3.057e-01   5.318 1.05e-07 ***
ethnicityMitte-eestlane  5.732e-03  2.471e-01   0.023 0.981495    
as.numeric(wage)         9.712e-02  3.544e-02   2.740 0.006141 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 770.46  on 795  degrees of freedom
Residual deviance: 710.26  on 788  degrees of freedom
AIC: 726.26

Number of Fisher Scoring iterations: 5

Kuidas visualiseerida interaktsioonide marginaalefektie (ehk tingimuslik tõenäosus)? Sellist võimalust pakub funktsioon interplot(...) paketist interplot.

library(interplot)
interplot(osalus_mudel, var1 = "age", var2 = "age") +
  theme_bw()

Rohkem võimalusi interploti kohta leiab siit: interplot: Plot the Effects of Variables in Interaction Terms

Et saada kätte lihtsamini tõlgendatavad šansside suhted, tuleb \(\beta\) väärtused eksponentida. \(\beta\) väärtused saame funktsiooniga coef() ning eksponendid exp(...).

(or = exp(coef(osalus_mudel)))
            (Intercept)                     age                I(age^2) 
              0.2150577               1.0319259               0.9999717 
            genderNaine    educationKeskharidus    educationKõrgharidus 
              0.9150965               2.3213569               5.0810969 
ethnicityMitte-eestlane        as.numeric(wage) 
              1.0057483               1.1019953

Šansside suhet ei saa siiski lineaarselt tõlgendada, sest konkreetne šansi suurus sõltub X väärtustest. Üks viis, kuidas kõikide X peale keskmist mõju suurust väljendada, on leida keskmine marginaalefekt. Selleks tuleb leida esmalt keskmine logaritm-suhte tihedus ning korrutada see koefitsientide suurustega. Kui see tundub praegu keeruline, siis antud hetkel ei olegi tarvis seda mõista. Huvi korral võib alustada tõenäosuse ja logaritm-suhte uurimist siit: Binary Outcome GLM Plots.

(ame = mean(dlogis(predict(osalus_mudel, type="link"))) * coef(osalus_mudel))
            (Intercept)                     age                I(age^2) 
          -2.163208e-01            4.423521e-03           -3.981973e-06 
            genderNaine    educationKeskharidus    educationKõrgharidus 
          -1.248868e-02            1.185380e-01            2.288028e-01 
ethnicityMitte-eestlane        as.numeric(wage) 
           8.067866e-04            1.367057e-02

Kui šansside suhet võis mõista “kui palju kasvab sündmuse šanss (\(\frac{y=1}{y=0}\)) X-i kasvades ühe võrra” (mitte segi ajada riskisuhtega (\(\frac{y=1}{y}\))!), siis keskmine marginaalefekt tähistab regressiooni joone keskmist kallakut.

Kuvame kõik leitud väärtused ühes tabelis kasutades selleks stargazerit.

stargazer(osalus_mudel,osalus_mudel,osalus_mudel,
          dep.var.labels = "Hääletas valimistel",
          column.labels = c("Beta koefitsiendid",
                            "Šansside suhe",
                            "Keskmine marginaalefekt"),
          dep.var.caption = "",
          coef = list(NULL, or, ame),
          type="html")
Hääletas valimistel
Beta koefitsiendid Šansside suhe Keskmine marginaalefekt
(1) (2) (3)
age 0.031 1.032*** 0.004
(0.032) (0.032) (0.032)
I(age2) -0.00003 1.000*** -0.00000
(0.0003) (0.0003) (0.0003)
genderNaine -0.089 0.915*** -0.012
(0.195) (0.195) (0.195)
educationKeskharidus 0.842*** 2.321*** 0.119
(0.244) (0.244) (0.244)
educationKõrgharidus 1.626*** 5.081*** 0.229
(0.306) (0.306) (0.306)
ethnicityMitte-eestlane 0.006 1.006*** 0.001
(0.247) (0.247) (0.247)
as.numeric(wage) 0.097*** 1.102*** 0.014
(0.035) (0.035) (0.035)
Constant -1.537* 0.215 -0.216
(0.788) (0.788) (0.788)
Observations 796 796 796
Log Likelihood -355.130 -355.130 -355.130
Akaike Inf. Crit. 726.259 726.259 726.259
Note: p<0.1; p<0.05; p<0.01

Arvutame viimaks McFaddeni pseudo R-ruudu suuruse. Arvestada tuleks, et erinevalt lineaarse regressiooni R-ruudust, on antud pseudo R-ruudu suurus alati madalam (0.2-0.4 on juba hea tulemus).

\(R^2=1-\frac{logLik_{mod}}{logLik_{null}}\)

1 - osalus_mudel$deviance / osalus_mudel$null.deviance
[1] 0.07813415

Meenutame, et andmestikus on hääletanuid oluliselt rohkem kui mittehääletanuid. Mudelit saaks seega oluliselt parandada tasakaalustatud andmestikuga. Kasutame selleks paketti ROSE ja funktsiooni ovun.sample(...).

library(ROSE)
over = ovun.sample(vote~.,data=RK_2015,method="over",N=1292)$data
under = ovun.sample(vote~.,data=RK_2015,method="under",N=300)$data
both = ovun.sample(vote~.,data=RK_2015,method="both",p=.5,N=nrow(RK_2015))$data

osalus_over = glm(f, data=over, family=binomial)
osalus_under = glm(f, data=under, family=binomial)
osalus_both = glm(f, data=both, family=binomial)

1 - osalus_over$deviance / osalus_over$null.deviance
[1] 0.08366703
1 - osalus_under$deviance / osalus_under$null.deviance
[1] 0.07123121
1 - osalus_both$deviance / osalus_both$null.deviance
[1] 0.1038426

Hea. Pseudo R-ruut tõusis tasakaalustatud andmetega mõne protsendi võrra.

8.3 Mudeli hindamine

8.3.1 Mudeli olulisus

  • \(\chi^2\) test

\(\chi^2\) test kontrollib, kas testidu mudel on oluliselt parem tühjast mudelist.

  • Hosmer & Lemeshow test

H&L test kontrollib, kas mudel sobitub halvasti. Hea on seega kõrge p-väärtus, mis näitab, et mudel sobitub hästi.

hoslem.test(osalus_mudel$y, fitted(osalus_mudel), g=10)

    Hosmer and Lemeshow goodness of fit (GOF) test

data:  osalus_mudel$y, fitted(osalus_mudel)
X-squared = 8.6897, df = 8, p-value = 0.3691

8.3.2 Klassifikatsioonitabel

logit_preds = ifelse(predict(osalus_mudel,type="response")>.5,"Hääletas","Ei hääletanud")
confusionMatrix(logit_preds,RK_2015$vote)
Confusion Matrix and Statistics

               Reference
Prediction      Ei hääletanud Hääletas
  Ei hääletanud             9        7
  Hääletas                141      639
                                          
               Accuracy : 0.8141          
                 95% CI : (0.7853, 0.8405)
    No Information Rate : 0.8116          
    P-Value [Acc > NIR] : 0.4496          
                                          
                  Kappa : 0.0748          
 Mcnemar's Test P-Value : <2e-16          
                                          
            Sensitivity : 0.06000         
            Specificity : 0.98916         
         Pos Pred Value : 0.56250         
         Neg Pred Value : 0.81923         
             Prevalence : 0.18844         
         Detection Rate : 0.01131         
   Detection Prevalence : 0.02010         
      Balanced Accuracy : 0.52458         
                                          
       'Positive' Class : Ei hääletanud

8.3.3 ROC-kõver

ROC_curve = function(probs,response,cut=.5){
  y_values = unique(response)
  cutoffs = sort(unique(as.numeric(probs)))
  tpr = sapply(cutoffs, function(cuts) sum((probs >= cuts) == T & response == y_values[1])/(sum(response == y_values[1])))
  fpr = sapply(cutoffs, function(cuts) sum((probs >= cuts) == T & response == y_values[2])/(sum(response == y_values[2])))
  stats = data.frame(cutoffs, tpr, fpr)
  
  r = rep(y_values[2],length(response))
  r[probs>cut] = y_values[1]
  cm = confusionMatrix(r,response)
  auc = cm$byClass["Balanced Accuracy"]
  
  stats %>%
    ggplot() +
    geom_line(aes(fpr,tpr)) +
    geom_ribbon(aes(x=fpr,ymax=tpr,ymin=fpr),alpha=.5,fill="steelblue") +
    geom_label(aes(.5,.5,label=paste("AUC:",round(auc,3))))+
    scale_x_continuous(breaks = seq(0,1,.1)) +
    scale_y_continuous(breaks = seq(0,1,.1)) +
    coord_fixed() +
    theme_bw()
}

ROC_curve(predict(osalus_mudel,type="response"), RK_2015$vote)